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ABSTRACT 


The primary goal of this work is to present for a control 
system a computer-aided-compensator design technique from a frequency 
domain point of view. Up to the present there have only been two simi- 
lar procedures developed, and they are somewhat limited. The thesis 
for developing the above said technique is to describe the open loop 
frequency response by n discrete frequency points which result in n 
functions of the compensator coefficients; several of these functions 
are chosen so that the system specifications are properly portrayed; 
then mathematical programming is used to improve all of these functions 
which have values below minimum standards. 

In order to do this several definitions in regard to measuring 
the performance of a system in the frequency domain are given, e.g., 
relative stability, relative attenuation, proper phasing, etc. Next, 
theorems which govern the number of compensator coefficients necessary 
to make improvements in a certain number of functions are proved. 

After this a mathematical programming tool for aiding in the solution 
of the problem is developed. This tool is called the constraint im- 
provement algorithm. Then for applying the constraint improvement 
algorithm generalized gradients for the constraints are derived. 

Finally, the necessary theory is incorporated in a computer 

program called CIP (Compensator Improvement Program). The practical 

usefulness of CIP is demonstrated by two large system examples. 

* * 
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' CHAPTER I 


INTRODUCTION 


Before the invention of the digital computer , elaborate and com- 
plicated numerical techniques for solving problems in science, mathe- 
matics, and engineering were only given secondary consideration*, As 
the refinement of the digital computer progressed, its comprehensive 
usefulness became more apparent. Today, the employment of the digital 
computer is found in almost every discipline of science and engineer- 
ing. 

Mathematical Programming 

One area in which the digital computer has been of tremendous aid 
is in the solution of mathematical programming problems. The general 
mathematical programming problem may be stated as: 

T 

determine the n component vector x - (x . , - x 2 , „ , 0 , x n ) so that 

T 

the maximum (or minimum) of f (x ) 

T 

subject to g^(x ){<_, = , >_} c^, 

i =1, 2,oco, m (1-1) 


is obtained. Each relation in (1-1) is assumed to be algebraic in 

T 

nature. The relation, f (x ), is called the , cost ; function, whose 
extremal with respect to the m constraints of the second relation is 
desired. If all the functions in (1-1) are linear and if the variables 
are not required to be Integral valued, then the above optimization 
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problem is said to be a continuous linear programming problem (LP). 

The solution of the continuous linear .programming problem may be 
accomplished with the aid of the simplex algorithm first introduced 
by George Dantzig. 1 Today the solution of continuous' linear pro- 
gramming problems is treated extensively in' many text books . 2 9 3 9 ** y 5 

On the other hand, if any of the algebraic functions in (1-1) are 
nonlinear, then the problem is called a nonlinear programming problem 
(NLP) 0 Up to now there has been no one algorithm developed that will 
solve all nonlinear programming problems. Generally the existence 
and uniqueness of a solution cannot even be assured without the cost 
function and the constraints possessing certain convexity and con- 
cavity properties o By placing various restrictions on the functions 
in (1-1) , there have been several algorithms developed for obtaining 
solutions o 6 In general, NLP algorithms are classed as either simplex 
in nature or as gradient in nature. 


Simplex Algorithms - Probably the- first NLP algorithms developed 
were the separable programming algorithms. Problems for which they 
are applicable are of the following form: 


n 


l gyCXjMl, ", >)C ± 


j"! 


X. 10 j = 1 , 


n 


maximize (or minimize) z = £ f ,(x.) 

j=i J 2 


( 1 - 2 ) 


Here dynamic programming is not considered as a NLP algorithm 
but is considered as another branch of mathematical programming . 
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In order to apply separable programming both the constraints and the 
cost functions must be separable into functions of' single variables. 
The mono-variable functions are then approximated over some finite 
interval by sequences of straight" lines. ‘ Then a simplex algorithm is 
used to solve the approximate problem. The separable programming 
algorithms differ in the way the approximations are made and in the 
type of simplex algorithm necessary to solve the problem, 6 

Another simplex NLP algorithm is the quadratic programming of 
Wolfe. 11 It was especially developed' to solve problems of the form: 


Ax = b 
x >_ 0 

T 

maximize (minimize) z = cx + x Dx 


(1-3) 


where A is an mx n matrix, c is an nx 1 matrix, and D is an n x n 
negative semidefinite matrix. In this case the constraints are linear 
and the cost function is quadratic and concave. The development of 
the algorithm for solving (1-3) depends heavily upon the Kuhn-Tucker 
conditions . 2 » 7 * 1 1 

Still another simplex type algorithm is the Hocking-Hartley con- 

1 2 

vex programming technique* It is used to solve general NLP programming 
problems with certain convexity and concavity conditions 0 This method 
is derived by . approximating the cost function and the constraints by 
an infinite number of supporting hyperplanes* Of course this produces 
a LP with an infinite number of rows » Then by using the duality prin- 
ciple of LP the problem is transformed into an infinite column problem 
which is amenable to solution by the simplex method. The convergence 
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properties of this algorithm are very reminiscent of the Newton-Raphson 
method for finding the roots of a polynomial, i. e. , whenever the 
algorithm converges, it usually converges’ very rapidly. 12 

Of course, there are many other simplex type NLP algorithms; in 
fact, there are several versions of those given above. However, for 
brevity only the more publicized algorithms and the' basic thoughts 
behind them have been mentioned here. 

Gradient Algorithms — In contrast’ to the- simplex type algorithms 
there exist the gradient algorithms. The premier algorithms of this 
type are the gradient projection method 1 3 > 11+ , the generalized reduced 
gradient method (GRG) 18 , and the sequentially unconstrained minimi- 
zation technique (SUMT) 1 5 » 16 » 17 

The basic idea of the gradient projection method is to start with 

a feasible solution and move in the direction of the gradient of the 

cost function (for maximization problems) until the solution is found 

* 

or until the violation of a constraint is attempted. If the viola- 
tion of a, constraint is attempted, a direction is determined so that 
an increase in the cost function results and' no' violation of the 
constraints occurs. If no direction can be determined then the 
solution has been found e 

In the case of linear constraints this simply requires projecting 
the gradient into the space defined by the intersection of all con- 
straints which are equalities at the point' under consideration. This 
is done by determining 

r -.Pd (1-4) 

A feasible solution is any point where no constraint is violated. 
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where r is the directional vector which points in the direction to 
move, d is the gradient of the cost function, and P is a projection 
matrix. The projection matrix is determined as 

P = I - Q(Q T Q) _1 Q T (1-5) 

where Q is a matrix whose columns are the gradients of the constraints 
which are strict equalities at the point of question. Of course if 
Q becomes square then P = 0. This does not indicate a solution but 
simply indicates that the feasible solution is located at a corner of 
the solution space. (For determining the projected gradient for this 
case, see Hadley [6], p. 167). 

Another so-called gradient NLP method is the GRG method mentioned 
previously. This technique is a natural extension of the reduced 
gradient method of Wolfe to include nonlinear constraints. The 
reduced gradient method was developed to determine relative extremals 
of 

maximize f(x) 

subject to Ax < b (1-6) 

■ x >_ 0 i = 1,2 , . . . ,n 

It is assumed that any n-row submatrix of A has rank n. Next, A is 
partitioned into an nx n submatrix C and a submatrix D, and b is 
similarly partitioned into c and d. Then slack variables y and z are 
added so that the constraints in (1-6) become 

Cx + y = c (1-7) 


Dx + z 


d . 


( 1 - 8 ) 
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All the constraints which are equalities are included in the C matrix. 
The variables of z are considered as dependent and' those of y 
independent. From (1-6), (1-7), and (1-8) it is easily seen that 



Ax 

■ = 

- C 1 Ay 

(1-9) 


Az 

= 

DC -1 Ay 

(1-10) 


y <*> 

= 

- Vf(x)C -1 

(1-11) 

where Ax 

and Az represent i 

the i 

changes in the x's and y's. 

V f (x) 

y 

is called 

the reduced gradient 

and 

Vf (x) is the gradient of 

the cost 

function. 

From (1-9), (1-10), 

and 

(1-11) a set of rules has 

been 


devised for determining the correct changes in the x T s and y T s so that 
an increase in f(x) is registered (For additional information see 
[30]). 

Somewhat different from the gradient projection and GRG methods 
is the SUMT. The problems amenable to this technique are those which 
can be cast into the following form: 

T 

minimize f (x ) 

T 

subject to g^x ) > 0 , i - 1, 2, . . , q 

h (x T ) = 0 , j « 1, 2, .... p . (1-12) 

In applying SUMT the above constrained minimization problem is trans- 
formed and solved as a sequence of unconstrained minimization problems 
which in the limit converges to a solution. This is done by forming 
from the above cost function and constraints a penalty function of 
the following form: 
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P(x T ,R) = f(x T )+R ^ — + -| f h 2 (x T ) (1-13) 

i=l g.tx 1 ) j=l J 

where R is a weighting constant greater than 0. For some initial 
value of R the unconstrained penalty "function, (1-13), is minimized 
by some unconstrained minimization technique. Then R is decreased by 
dividing it by some number greater than 1 and the process is repeated. 
As R -*■ 0 the unconstrained solution approaches a constrained solution. 
The physical effects of the two latter terms in (1-13) is to penalize 
a trial solution for getting too close to the boundary of the feasible 
region. 

There has been no attempt here to be all inclusive with respect 
to gradient algorithms. There are several other- gradient algorithms 
that have been developed. However, the ones mentioned above are con- 
sidered by many as the most prominent and useful methods . today . 

Mathematical Programming in the Design of Control Systems 

Over the past ten years there has been a great thrust to use 
mathematical programming in the design of control systems. The major 
effort has been in the solution of optimal control problems , and the 
results in this area have been very fruitful* — not only in the appli- 
cation of mathematical programming but also in theoretical develop- 
ments, In fact, it has been shown that the Kuhn-Tucker necessary 
conditions of mathematical programming and the maximum principle of 
optimal control can be derived from the same set of general optimiza- 
tion theorems 19 * 20 * 21 ’ 22 . As can be seen from the lengthy reference 
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list by Tabak 23 , much of the work has been directed toward the appli- 
cations of linear and quadratic programming. Recently, uses of the 
SUMT and the GRG algorithms in the solution of optimal control 
problems have been made. 2i+ > 25 > 26 

On the other hand, the use of mathematical programming in the 
classical design of control systems has been meager-^particularly in 
the design of compensators from a frequency domain point of view. 

This is very unfortunate because most practical system designs even 
today are still by classical frequency domain approaches. Further- 
more, these approaches are more artful than analytical. The few 
techniques which have been developed can be classified as modern con- 
trol oriented or strictly classical control oriented. This classifi- 
cation results from the choices of the performance indices. Those 
methods in which system specifications are submerged in a cost 
functional are labeled as modern control approaches, while those 
methods , which represent the system performance by classical standards 
such as gain margins, phase margins, bandwidth, etc * , are termed as 
classical approaches . 

One of the first successful computerized compensator algorithms 
was developed by Coffey. 27 In his paper - consideration is given to 
a system similar to that shown in Figure 1. In this figure j parame- 
ters of the system are sensed; each parameter is operated on by some 
compensation device; the results of these are summed and fed back. 
Figure 1 is considered typical of large aircraft or space vehicles. 
Each compensator is assumed in the following form: 
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M N 

G (s) = [ f a fb .s 1 ' 1 ] (1-14) 

e i=l el i=2 ei 

where s is a complex variable and - 1 and N^ - 1 are the nume- 

th 

rator and the denominatororders'; respectively , of the e* compensator. 
The goal is to select the compensator coefficients so that the com- 
pensated open loop frequency response' is' a weighted least squares fit 
to a desired open loop frequency response (Of course, the open loop 
frequency response is obtained by calculating C (jo))/R(jo)) when the 
feedback loop is broken at a). 

The weighted least-squares fit is obtained by minimizing the 
following cost function: 

J “ I (y* - y*) T WW T (y - y) | (1-15) 

where y is a vector of the desired frequency response points, y is a 
vector of frequency response points, W is a diagonal weighting matrix, 
the asterisk (*) denotes conjugate, and the super T denotes transpose. 
For minimizing the cost function, J, with respect to the compensator 
coefficients a gradient search algorithm is chosen; and, then, the 
necessary gradient vector is calculated. 

Next, geometrical properties of the cost ; function are considered. 
It is demonstrated that even for relatively simple systems the cost 
function is geometrically complicated. From this it is seen that the 
cost functions can have relative extremals and unbounded solutions. 
Furthermore, the design of unstable compensators is possible. Even 
with the possibility of these difficulties, it is demonstrated that 
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this procedure can be utilized to- design practical compensators . This 
is done by applying the technique . to* a sixth order ballistic missile 
example. For this system two compensators are designed—a pure gain 
and a fourth order over a sixth' order. The pure gain compensator 
approximated the desired frequency response for low frequencies but 
was completely unsatisfactory for'higher frequencies. In fact, for 
this compensator the closed loop system" is unstable. On the other 
hand the higher order compensator exhibited very good properties when 
compared to the desired frequency response. 

Coffey indicates that in some instances a judicious choice of 
the elements of the weighting matrix, W, is required before an 
acceptable design can be achieved. Thus, a computer program of this 
algorithm might require several runs^-while juggling these elements 
between runs — before the proper values are conceived. Even with this 
disadvantage the algorithm is definitely superior to classical means. 

Another technique for computerized design of compensators for 
control systems has been presented by Page and Stear e 28 » 29 The thesis 
of this procedure is to vary the compensator coefficients until 
certain chosen frequency response specifications are satisfied. The 
procedure for attempting to do this is 

N , 

minimize F = £ K. (1 - S^/S. ) (1-16) 

i=i 1 11 

a 

where N is the number of specif ications considered , is the speci- 
fication as a function of the compensator coefficients, is the 

desired specification, and is a weighting constant. The constant 

a d 

is chosen as positive, in general one, for 


and as zero 
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a H 

for . This results in a satisfied specification being 

neglected. The goal is to drive F to zero. The reason for the choice 
of the above criterion function (1-16) is to try to place the most 
emphasis on the specifications which' have the greatest violations. 

In order to illustrate the given procedure Stear and Page pre- 
sent an example of the design of an autopilot for an aircraft. In 
accomplishing this design four unconstrained optimization procedures 
are used. Three are local search procedures, and one is a global 
search technique. As in the case of Coffey's cost function it Is 
discovered that even for simple compensators the specification 
function (1-16) has relative extremals. From this it is deduced that 
the global search procedure is more applicable than the local search 
techniques if the starting compensator is strictly arbitrary. However, 
if a priori knowledge is used in picking the initial compensator this 
deduction is not necessarily true. 

Pitfalls of Previous Works on ' Computerized ' Compensator Design 
Procedures 

The two previously mentioned works on computerized compensator 
design procedures suffer from several drawbacks. First the procedure 
presented by Coffey is basically a frequency response shaping technique. 
In the design of compensators for most control systems, this is too 
rigorous; i. e. , this requires the compensator to satisfy more con- 
straints than are necessary. Thus, the probability of all system 
specifications being satisfied is less. Another interesting fact is 
that in many instances the frequency responses of control systems 
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are not required to match a desired frequency response—frequency to 
frequency — but are desired to have some general shape which can be 
translated with respect to frequency. Even more conceivable is the 
desirability to have several bands of the frequency response to be 
various distances from the -1 + jO point of the GH(jw)-plane and to 
have other bands of the frequency response constrained to be greater 
than or less than limitations with respect to the origin of the 
GH(ja)) -plane. Constraints such as these are not as strenuous as those 
requiring the frequency response to fit closely to some desired fre- 
quency response. 

A pitfall which is common to both the Coffey method and the 
Stear and Page method is the necessity of choosing some constants — in 
particular, the elements of the diagonal matrix, W, and the K_ f s. It 
is obvious that in many situations a judicious choice of these must 
be made before any useful results will emerge. It was suggested by 
these authors that computer programs containing the algorithms may 
require several runs with various values of these constants before an 
acceptable design is achieved. However , this involves trial and 
error which was one of the justifications for going to a computerized 
procedure. 

Another drawback of the two algorithms presented is that some 
specifications may become worse while others become better . This 
immediately poses some serious questions, such as, what is a reason- 
able trade-off and where does it exist? If minimum standards of 
system performance have been set, it is very probable that nothing 
short of these are acceptable. In this case there is no trade-off. 



14 


On the other hand, it may be viewed' that - in' practical designs it is 
not unusual to accept performances a little less than that desired. 

In instances such as this, performance tolerances must be set. 

Another shortcoming of the two methods is their failure to 
include inherent devices for maintaining compensator stability. If 
the designed compensator is unstable, then the stability criterion of 
the system changes completely. The result might be system instability 
which removes the compensator from the realm of a practical design. 

What is needed is an algorithm which tends to improve system specifi- 
cations at every iteration. Of course this might require the allowance 
of only incremental changes in the compensator coefficients. 

Another pitfall of the two previously mentioned works is the lack 
of any theoretical inclusions on compensator limitations . That is, 
none of the authors presented any theoretical developments showing 
what could be. expected from their algorithms for a certain compensator 
order in a particular system. Thus, initially there is no way to know 
what minimum amount of compensation is necessary. In addition, these 
works presented no theory which' indicates that the algorithms will 
produce a final compensator that’ is any better than the initial com- 
pensator. 

In essence, the techniques of Coffey and Stear and Page are 
"firsts" in the use of the computer for compensator designs, but they 
are somewhat limited. They do not present universal solutions in 
regard to computerized-compensation. It is the purpose of this dis- 
sertation to present the theory and a method of computer-aided 
compensator design that does not have the drawbacks of the previously 
presented techniques and is thus more. universal. 



CHAPTER II 


FREQUENCY RESPONSE CONSIDERATIONS IN THE 
DESIGN OF' A CONTROL SYSTEM 

Before the design of a system can be accomplished, the limitations 
or constraints and the desired performance of the system must be 
established. The measurement of the performance of the system is 
determined by comparing that obtained to that desired. Because of the 
limitations, in many instances, the desired performance cannot be 
achieved. In designing compensators for practical control systems 
there are, in general, two types of performance indices- — time domain 
indices and frequency domain indices. Although it is quite obvious 
that these are related, no analytical means, up to now, have been 
devised for defining this relation except for the simplest control 
systems — less than third order. In practical designs the main limi- 
tations are system stability, nonlinearity, time variance, and 
sensitivity. Today many systems are designed by' using linearized 
frozen time models and applying frequency domain concepts. 

Concept of Relative Stability 

In most practical systems stability' is a' major constraint. In 
fact, inmost system designs a specified degree of stability is 
required. A specific degree of relative stability is required because 
of inaccuracies in the model of the system or in order to deter insta- 
bility if future parameter variations in the system plant result. 



16 


Sometimes a certain amount of relative stability is desired to keep 
the system from resonating unnecessarily., 

In the past the degree of relative stability has been denoted 
by the classical gain (GM) and phase' margins (PM) « However, in some 
instances these can be very misleading. For example, consider the 
hypothetical s-plane frequency response' shown in Figure 2 which 
possesses acceptable classical stability margins (GM >_ 2.0, PM >_ 30°) 
but which comes within some small distance of the -1 + jO point. Such 
a condition could represent a system which was, very close to insta- 
bility. A better measurement of relative stability is defined as 
follows : 

A stability margin is defined as the magnitude of the 1 + GH(joj) 
frequency response ;at one of its minima relative to the origin 
of the 1 + GH(ja)) plane. 

It is deemed by this author that by' measuring stability in this 
fashion, a measure of the true relative stability of a system is 
achieved. Next, a system is said to be relative stable if the fre- 
quency response does not cross ~ a designated closed contour located 
around the -1 + jO point. This closed contour around the -1 + jO 
point is called the margin of stability limit. 7 - The shape and the 
size of this contour depend upon system specifications. Furthermore, 
there is nothing wrong with making the size and shape; of the contour 
frequency dependent. (In doing this the designer would be indicating 
that the frequency response is to be shaped to some extent.) 

Relative Attenuation Concept 


Although relative stability plans a major role in compensator 
determination, there are several other factors which are considered. 


17 



Figure 2. A Hypothetical GH(joj) Frequency Response 


One of these is the attenuation of certain frequency bands. The reason 
for frequency band attenuation is to discourage the control system 
from resonating at some natural frequency of the system. Of course if 
the system is linear and time-invariant this is not necessary. Un- 
fortunately, many practical systems do not fit* into the linear , time- 
invariant category. 

Frequency band attenuation may be* treated by requiring that all 
frequency points that are to be attenuated fall within a chosen con- 
tour around the origin in the GH(joj) plane. This contour is called 
the margin of attenuation limit . It then follows that: 

An attenuation margin is the magnitude of the GH(jw) 
frequency response at one of its* maxima with respect to 
the origin of the GH(ja)) plane . 7 

Other Frequency Response Concepts 

Relative stability and attenuation are considered as the most 
important frequency response design criteria. However, they do not 
yield acceptable designs in all instances. Sometimes it is necessary 
to employ proper phasing of certain frequencies. This is usually 
employed when it becomes difficult to determine a compensator to 
attenuate certain natural frequencies' of the system and in, addition 
to satisfy other system requirements. The general idea is to 
determine a compensator so that these frequencies are phased toward 
the right half of the GH(jca) plane. This results in these frequencies 
being attenuated in the closed loop system. 

In some cases it is even* necessary to place special emphasis on 
certain points of the frequency response. In most instances these 
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points are closely related to' dynamical' responses of the control 
system. Examples of dynamical' responses' considered for a space craft 
are wind response and "engine-out” response. In order that these 
responses possess acceptable characteristics it is usually necessary 
to require certain frequency response' points- to be placed in certain 
regions of the GH(j to) plane. 

I 

Still another frequency response design- concept is bandwidth • 
However, this can be handled by either the stability margin or the 
attenuation margins . For example, the maximum open loop bandwidth 
can be achieved by requiring a certain frequency and all frequencies 
above it to have a certain margin of attenuation limit. Similarly, 
closed loop bandwidth could be controlled by a combination of these. 

Problem Formulation 

Assuming that the desired frequency response characteristics have 
been determined so that if they are achieved the performance of the 
system will be acceptable, it must be decided how to determine a 
compensator for achieving these. The classical means of doing this 
is by trial and error; however, a more efficient method would be an 
iterative method that makes improvements upon the system* s frequency 
response from iteration to ; iteration or indicates that no further 
improvement could be made. In fact, if a total of n critical frequency 
points have been chosen, then the problem may be formulated as the 
following nonlinear programming problem: 

T 

Determine a vector x such that 
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g ± (x T , (op . <. b i 

: ", ' =, • >} d i 

i = 1, . . . , n (2-1) 

T 

In (2-1) x is a vector of the compensator coefficients; g^ is a 
function of the i^ frequency, o>^, and the* compensator coefficients. 

The functions, g^, i = 1, ..., n~, are' chosen so as to represent the 
frequency response limitations and constraints which have been imposed. 
For example, g^ could be representative of a. stability margin or an 
attenuation margin. The second relation in (2-1) takes into account 
any constraints that might be placed on the compensator coef f icients. 

It may be necessary to constrain some of the coefficients if it is 
desired to keep the d. c 6 gain, G(jO), of the system constant or above 
or below a certain level. Also, it may be necessary to constrain 
certain compensator coefficients to insure the stability of the 
compensator or to take into account realizability conditions. 

The above formulated nonlinear programming problem differs from 
the classical nonlinear programming problem in the respect that it is 
strictly a constraint problem. 6 * There is no cost function to maximize 
or minimize. However, this does not simplify matter s 0 In fact, the 
above problem can be thought of as a normal nonlinear programming 
problem in which it is desired to find a solution which obtains a 
certain objective function value. In this case the objective function 
just becomes a constraint. If the objective function is added to the 
constraint list, then the result is a strict constraint problem as 
given above. The desired solution to this problem is a feasible 
solution which may or may not existo 



CHAPTER III 


• COMPENSATOR LIMITATIONS 

At any iteration in' solving the' problem mentioned in Chapter II, 
there will result conditions of the form of (2^1) to be improvedc 
(The number n can change from one iteration to another since the 
frequency response changes with respect to the compensator . ) The 
general idea is to change the compensator coefficients so that each 
constraint comes closer to being satisfied. The question then is, 

how many compensator coefficients are required to insure that some 

i 

improvement on each constraint at a certain iteration can be made? 
This question is answered by the following definitions. and theorems. 

Definition 1 

An optimal direction in the GH(joj) plane is any chosen 
direction in which it is desired to perturb a point on 
the frequency response. 

Optimal directions are illustrated in Figure 3 at points A, B, and 
C. The number of compensator coefficients sufficient to perturb n 
polar frequency response points in n optimal directions is given by 
the following theorem: 

Theorem 1_ 

A sufficient condition to' perturb n points on a polar frequency 
response curve in n optimal directions with a realizable , compensator 
is that there be at least 2n independent compensator coefficients 


which are available to be varied. 



GH(Jco)-PLANE 



Figure 3. A GH(joo) Frequency Response for Illustrating 
Optimal and Sub-optimal Directions 


23 


Proof : Let the open loop frequency response be denoted by 

T T 

G (jw, x ) where x is an m dimensional vector of the functionally 
o 

independent compensator coefficients. Also, let the optimal 

* 

direction at a frequency ou^ be denoted by d^ . Suppose there are n 
points on the frequency response' which are to be moved in the n 

chosen directions, respectively. The change of the open loop transfer 

th th 

function at the 1c frequency with" respect to the i" compensator 

coefficient is of the form 


3G o (jco k , x 1 ) 


3x , 


C ki + je ki 


(3-1) 


where and e^ are real constants. There are, for a particular 
frequency, m such partials as (3-1) and, if they were included as the 
components of a single vector, the result would be the complex 
gradient. It is well known that this points in the direction of the 
most rapid change. However, this is not the desired direction of 
movement. Essentially what is needed is a directional vector [w] in 
complex m-space whose dot product with the m dimensional gradient 

*J5* 

vector [c^. + je^] will yield the desired directional derivative d^_ , 
or in equation form (See [32]) 


d k = [c k + je k ] [w] 


(3-2) 


It should be obvious that the components of [w] are proportional to 
the amount that each compensator coefficient must be varied in order 

k 

that movement in the d^ direction can be accomplished. Thus if the 
compensator is to be realizable, [w] must be a real vector. 
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Letting 



(3-3) 


then (3-2) can be written by the following two real equations: 


a k * = [c k J T [wj (3-4a) 

and 

b k* “ [e k ]T [w] • (3 " 4b) 


Hence, for n points on the frequency response to be moved in n 
optimal directions there result 2n equations or 

aj = [c 1 ] X [w] 

•k T 

b 2 = [c 2 r [w] 


* 

a = 

[c] T 

[w] 

n 

n 

* 


[w] 

b i = 

te i ] 

* 

b 2 = 

[e 2 ] T 

[w] 

* 

b 

n 


[w] . 


In matrix notation (3-5) becomes 


[w] 


(3-6) 


where the dimensions of 


* 

a 

• • • • 

„ * 


c 

• • • < 


T 

L e 1 J 


, and [w] are 


respectively 2n x 1, 2n x m, and m x 1. If 2n > m there will result 
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more equations than unknowns and possibly an incompatibility. 31 Hence 
there may not exist a vector [w] such that all equations can be satis- 
fied. This says there are not enough compensator coefficients avail- 
able. On the other hand if 2n <_ m, there either results less equations 
than unknowns or the same equations as unknowns. For the first case 
there will exist an infinite number of vector [w] f s and an infinite 
number of solutions to the equations. This indicates an excessive 
number of compensator coefficients. In the - second case there will be 
a unique [w] and, thereby, a unique solution for the equations. This 
means that the exact number of compensator coefficients necessary is 
being employed. 8 

The preceding proof has shown the sufficiency condition for mov*- 
ing the frequency response in n optimal directions. Suppose, however, 
that it is desirable to use a compensator with' a fewer number of 
coefficients than those needed to move' in the optimal directions. 
Consider the following definition: 

Definition 2 . 

A sub-optimal direction is any' direction within tt/2 

radians of an optimal direction. 

An optimal direction is just a two-space y ‘ vector ; then, a sub-optimal 
direction is any two-spacevector which has a positive dot product 
with an optimal direction. Thus, a sub-optimal direction is any 
vector which falls within a certain open half space, e.g. , a sub- 
optimal direction to B in Figure 3 is any vector which points to the 
left of the line passing through B. 

If the optimal and sub-optimal directions for are respectively 
represented in 2-space by the following vectors : 
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(3-7) 


and 


d k = (a k * V ’ (3_8) 


then the sub-optimal direction would be any direction such that the 
dot product 

d fc • d fc * > 0 (3-9) 

or 

a k a k* + b k\* > G • (3_10) 


Then the question is, how many compensator' coefficients are necessary 
in order to assure that movement in some sub^optimal direction can be 
achieved? The answer to this is stated and proved in the supervening 
theorem. 


Theorem 2 


In order to be assured of perturbing n points of an open loop 
frequency response in n sub-optimal directions, by varying the compen- 
sator coefficients, it is necessary- that there be n independent 
compensator coefficients available for variance, 

th 

Proof : The components of the k sub-optimal vector direction in 

terms of the real and imaginary parts of the partials at the 
frequency are given by 


m 



c 1 . w, 
ki i 


(3-11) 
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m 

b k ■ JjVi 


(3-12) 


where c^ and e^, respectively, are the real and imaginary parts 

(evaluated at w^) of the partial of the open loop transfer function 

til til 

with respect to the i compensator coefficient, and w^ is the i 

unknown constant which is to be determined so that (a^, b^) points in 

a sub-optimal direction. Substituting (3-11) and (3-12) into (3-10) 

results in 


I 


i=l 


'ki 


w. 



+ 


m 


l 

i=l 


e ki W i \ 


> 0 


(3-13) 


or 


m 


(c ki a k + v b k > ”i * 0 • 


(3-14) 


Remembering that there are n frequency points, n inequalities 
like (3-14) will result. Hence the following matrix inequality can 
be obtained: 

T 1 fp ^ 

[c a 4- e b ] [w] > 0 . (3-15) 

T * T * t 

The dimension of [c a + e b ] is n x m. In order to be assured that 
all n inequalities can be satisfied, it is necessary that there be at 
least the same number of unknowns as inequalities. Hence, this says 
there must be at least n independent* compensator coefficients in order 
to be assured that n frequency points can be perturbed in the sub- 
optimal directions. 8 

The above, two theorems place limitations on the overall compen- 
sator order. Thus for any algorithm to be assured of being able to 



make the changes given in the theorems, the theorem must be sat- 
isfied. 



CHAPTER IV 


CONSTRAINT IMPROVEMENT ALGORITHM 

It is very desirable to have an algorithm which starts with 
some initial compensator and, then, in an iterative fashion produces 
an improved frequency response. This statement immediately suggests 
the question — what is an improved frequency response? This is 
answered by the following two definitions. 

Definition 3_ 

A total improved frequency response (TIFR) in an iterative 
scheme is one whose unsatisfied constraint values at a 
certain iteration are better than they were at the last 
iteration. 

Definition 4_ 

A sum improved frequency response (SIFR) in an iterative 
scheme is one whose sum of the differences in the unsatis- 
fied constraint values and their desired values is a positive 
value from one iteration to the next.* 

It is obvious that an algorithm which is capable of producing a TIFR 
is also capable of producing a SIFR; however, this statement is not 
reversible. A TIFR algorithm requires every constraint which is 
unsatisfied to be improved or bettered at every iteration* while a 
SIFR algorithm only necessitates a sum improvement, i. e. , the sum 


It is assumed, here, that all constraints in (2-1) have been 
represented in the <_ form by multiplying >_ constraints by -1 and 
changing = constraints to two inequality constraints (See Hadley 
[6]). No generality is lost by doing this. 
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increase must be better than the sum decrease. The goal is then to 
derive an algorithm which is compatible to both TIFR and SIFR. 

Thus, an algorithm is needed for solving a nonlinear program- 
ming problem of the following form: 

T 

Determine the vector x such that 

T 

g.(x ) >. b ± i = 1, , m . (4-1) 

Again this is strictly a constraint problem. If this problem has a 
solution, then it is a point in a solution space (Theoretically the 
solution space could be a single point). The functions in (4-1) are 
not assumed either concave or convex. What is desired is an iterative 
algorithm which, when started at some initial guess at the solution, 
will at each iteration produce an improved solution from the solution 
at the last iteration or will indicate that no further improvement 
can be made. An improved solution is defined as one which brings the 
constraints closer to being satisfied. 

Constraint Improvement Algorithm Derivation 

T 

Suppose that some initial starting point, , has been chosen* 
Of the m constraints, let n be the number not satisfied by this point. 
The constraints not satisfied are defined as the active constraints 3 
and those satisfied are called the inactive constraints. Let J 
contain the index numbers of the active constraints, i. e. , 

J = (ki, k 2 , ...» k^}. Essentially what is desired is a directional 
vector, D, by which the vector x can be changed, and it will be possi- 
ble to get an improved solution. This vector can be calculated as 


D 


VSkj + a 2 V 8 k 2 + • • • + a n V gkn . 


(4-2) 
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In (4-2) 


k 


1 * 



k n e J 


9 


denotes the gradient of the constraint corresponding to the k^ 
index evaluated at , and {a^} is a set of constants that are to be 
determined. An improved solution can be assured if the a's are 
determined so that 



D • Vg k ^ > 0 i — 1 n 


(4-3) 


T 

In other words the maximum rate of increase of g, at x, is in the 



in the direction of any vector which has a positive component in the 
direction of the gradient. In fact, suppose that a value for each of 
the dot products in (4-3) is chosen. Then (4-3) becomes 


D * v gk 2 = C 1 
D • Vg k = c 2 

D ‘ V Sk n = c n 
T 

where the vector c = > c 2 * * • • * c n ) contains the chosen dot 

product resultants. Substituting (4-2) into (4-4) results in the 
following set of linear equations, 

< 7 8 kl * VSk^ 3 ! +-< v 8k 1 * V Sk 2 )a 2 + ••• + ( v g kl * V Sk n K " C 1 

( v 8k 2 * v 8kj> a l + (V 8k 2 * v Sk 2 > a 2 + ••• + (Vgk 2 * ? gkn )a n - c 2 


(V gkn • Vg ki ) ai + (Vg kn • Vg k2 )a 2 + ... +.(Vg kn - Vgj^Jan = c n . (4-5) 
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Using matrix notation (4-5) becomes 

[VG T VG]a = c , (4-6) 

T 

in which a = [a^ a 2 ... ajj] and VG is a matrix whose columns 

are composed of the gradients of the active constraints (The matrix 
T 

[VG VG] is the Gramian matrix of the gradient vectors under con- 
sideration — see Hildebrand [31].). 

If the gradient vectors are linearly independent then 

T “1 

a = [VG VG] c . (4-7) 

Hence, this will yield a's for a desired dot product between the 
directional vector D and each gradient of the active constraints. 9 
By moving in the direction of D then it is possible to improve the 
present solution.* 


Algorithm Summation 


Using the derivation and the preceding terminology, the con- 
straint improvement algorithm may be summarized as follows: 

x k+i ’ \ + h[VG)a 

T T tli tli 

in which and ara the solution points at the (k + 1) and k 


In the above derivation the gradients were used. However, 
vectors in the directions of the gradients will suffice. In fact, it 
has been found in practice that unit vectors in the directions of the 
gradients are more suitable when the gradient magnitudes become 
disproportioned. The main advantage is a greater convergence 
rate. 
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iterations respectively, [VG] is a matrix whose . columns are composed 

T 

of the gradients of the active constraints evaluated at x^ , 

T — 1 

a = [VG VG] c 

where c is a column matrix of positive constants, and h is a positive 
constant. 

The choice of h (the step size constant) determines how much or 
whether any improvement in the constraints is made. In a compensator 
design program h also determines whether the program is a TIFR or SIFR 
algorithm. As a general rule small positive values of h produce a 
TIFR and larger values of h produce a SIFR. Of course there is a max- 
imum limit on h for producing a SIFR, i. e. , values of h above the 
maximum do not produce either a TIFR or a SIFR. On the other hand, 
negative values of h are out of the question since they tend to 
decrease the constraints — making them even worse. 

In addition to choosing h, a choice of the components of the c 
vector must be made. As has been pointed out previously, the com- 
ponents of c are the dot products of the directional vector, D, and 
the gradients of the active constraints. Thus by properly choosing the 
c T s the amount of increase in some of the constraints can be, to some 
extent, controlled. In other words by judicious choice of the c*s some 
constraints can be weighted more heavily than others. However, the 
actual amount of change in a constraint is related to h and the con- 
straints* partial derivatives. In practice it has been found that 
when using unit vectors in the directions of the gradients of the con- 
straints a good choice of the elements of the c vector is l*s. This 
choice gives the best convergence rate. 
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On the other hand, there is nothing wrong with making the c T s 
dependent upon the constraint values, e. g., by letting a c decrease 
as its constraint comes closer to being satisfied. However, as a c 
approaches zero the algorithm would tend to determine a direction that 
was parallel to the boundary of the feasible region. Hence, the proba- 
bility of the constraint corresponding to this c becoming inactive 
decreases. Nevertheless, it has been discovered that in many instances 
that by holding the c T s at respectable positive levels many of the 
constraints are driven to inactivity and they do not return to activity 
again. In this case the order of the matrix whose inverse is required 
can be reduced, whereas, if all constraints always linger in activity 
the order can increase if other constraints become active on higher 
iterations. 

Algorithm Limitations and Termination 

Next, attention is focused on algorithm termination. There are 
three conditions in which the algorithm will terminate. These are 

1. All constraints are inactive. 

2. One of the gradients of one of the constraints becomes zero. 

3. The gradients of the active constraints become linearly 
dependent. 

The first of these simply indicates that' a solution has been obtained. 
The second and third represent relative extremal solutions . In fact, 
the second one shows that the solution point is a relative extremal of 
one of the constraints. On the contrary, the third termination con- 
dition indicates that at least one of the constraint gradients is a 
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linear combination of the others' gradients or there are more active 
constraints than there are variables (This could represent an incom- 
patibility condition.). Whenever 2 or 3 occurs either the solution 
obtained will have to be accepted or a new starting point will have 
to be chosen and the algorithm reinitiated. 



CHAPTER V 


GENERALIZED PARTIAL CALCULATIONS 

In essence, the goal of the designer is to pull and push various 
points on the frequency response until system specifications have been 
met or until no further improvements can be accomplished by the present 
compensator. In general, this can be accomplished by pushing and 
pulling the various points with respect to other points in the complex 
GH(jw) plane. For example, relative stability can be obtained by push- 
ing the points , of the stability margins away from the' -1 + jO point. 

On the other hand, the attenuation margins can be improved by pulling 
these points toward the origin. Similarly, proper phasing could be 
achieved by attempting to pull or push these points with respect to 
real axis points. Of course, in some specialized cases it may even be 
advantageous to pull or push a point' with respect to more than one 
point. Regardless of whether a point' is to be pushed or pulled it is 
necessary to know how these points change with respect to other points 
in the GH(jw) plane. This is especially true if the algorithm in 
Chapter IV is to be used in perturbing these points. 

A point can be pushed or pulled with respect to another point, 

- K, in the complex GH(joo) plane by varying the distance squared, 
d(w), between the point and ■- K. In order to determine how this dis- 
tance changes with respect to the compensator coefficients, con- 
sideration is given to the general feedback control system shown in 
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Figure 1. The open loop frequency response of this system is deter- 
mined by breaking the feedback loop at a and then calculating 

GH(jw) = C(jw)/R(jto) . (5-1) 

Furthermore, to generalize even further- in Figure 1 each channel's 

compensator is assumed to be made up of a product of sub-compensators, 
tVl 

i.e., the k channel's compensator is given as 

n k 

G,(s) = n G, . (s) , (5-2) 

i=l 

til ^ 

where n^ is the number of sub-compensators in the k channel. The 

th 

uncompensated open loop state frequency response of the k channel 
with all channels opened is defined as 


(j^) = a k (m) + jb k (u>) 


(5-3) 


where a^ is the real part and b^. is the imaginary part. 

From the above equations and statements it then follows that 


d (to) = 


J k 

K + l {[a k (u) + jb k (m)][ H G ki (jm)]} 
k=l i=l 


(5-4) 


By assuming each sub-compensator to be a general rational function of 
the following form 


G 

qp 


(s) 


n , n-1 . , 

xs 4-x .s + . . . + x^ 
n n-1 0 

m . m-1 . . * 

y m s +y m-l s + *-' +y 0 


(5-5) 


This is called the factored form of a , compensator . 
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it becomes necessary to derive only how d(w) changes with respect to 
the coefficients of this general compensator, because the change in d(m) 
with respect to any compensators’ coefficients will assume the same 
general form, only differing by the orders, n and m, and the numerical 
values of the x’s and y’s. Since G^(joj) is completely independent of 
all the other compensators, then it may be isolated from the others in 
(5-4), This is easily done by letting 

J n k 

A + jB = K + l {[a (o>) + jb, (u)][ n G, (ju) ] } (5-6) 

k=l i=l 

k^q 

and n 

q 

C + jd = [a (u) +-jb (w)] n G.(jw) . (5-7) 

4 4 i=l 4 1 

i^P 

Using (5-6) and (5-7), (5-4) is rewritten as 


d(u>) 


A + jB + (c + jd)G 

qp 



(5-8) 


Substituting (5-5) into (5-8) and carrying out the necessary manip- 
ulations (5-8) evolves as 


n 


( i-o Cl>tl + A j=0 E2j!,2j " B j=o E2j+1 yy+1> + 


( I D X + A I E y + B j t , ) ! 

d(.) - -i=2-j J -°- (5-9) 

( jJ 0 E 2j 7 2j ) + ( J 0 E 2)+l y 2j+l ) 
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where in (5-9) k = m/2 and p = m/2 - 1 
if m is even 

or k = (m-l )/2 and p = (m-l )/2 

if m is odd; 

the C's, D's, and E's are defined by the following sets: 


{ C o , C j j C 2 ) C 3 , C 4 , C 5 j * * • } ( c j doi ) — cgj 2 , dw^ ^ cto j — dco , •••} 

{Dq, D 1 , D 2 , D 3 , D 4 , D 5 , ...} = {d, cto, -dw 2 , -cu 3 , doj 4 , cgj 5 , ...} 

{Eq, El, E 2 , E 3 , E 4 , E 5 , ...} = {1, GJ, -GJ 2 , -GJ 3 , GJ 4 , GJ 5 , ...}. 

( 5 - 10 a,b,c) 


Next, letting 


FN1 = 


L C i X i + A Jq E 2 j y 2 j B E 2 j +1 y 2 j+l 


n 

l 

i =0 


FN2 = 


Jo Vi + A Jo E y +1 y « +1 + B Jo E « y « 


FD1 = 


Jo y « 


E 2 j +1 y 2 j+l 

2 2 

FD = (FDl) + (FD2) 

2 2 

FN = (FN1) + (FN2) 


then 


9d (gi) 
8 x 

q 


2[FD (A • FN1 + B • FN2) - FN • FDl] E q 
(FD ) 2 


(5-11) 

(5-12) 

(5-13) 

(5-14) 

(5-15) 

(5-16) 


(5-17) 
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for q even or 


3d (to) 2 [FD (-B » FN1 + A. • FN2) - FN. ♦ FD2]E q 
3x q (FD) 2 


(5-18) 


for q odd and 


3d(o>) 2 [FN1 * -C q + FN2 - D q ] 

3y q FD 


(5-19) 


for q even or odd. 

By programming equations (5-4), (5-6), (5-10a,b,c), and (5-11) - 

(5-19) on the digital computer the partials of d(aj) with respect to 

the coefficients G (s) can be obtained. 9510 
qp 

The above derivation provides the key for determining how any 
sub-compensator affects d(oj) in a first order sense. With a complete 
comprehension of this derivation it becomes clearly apparent how to 
proceed either from channel to channel or from sub-compensator to sub- 
compensator in order to determine the necessary partial derivatives 
for a particular frequency point. Of course, this process must be 
completely repeated for each individual frequency point. Once the 
gradient vectors of each chosen frequency point are determined, then 
the calculation of the directional vector is accomplished as described 


in Chapter IV. 



CHAPTER VI 


COMPENSATOR IMPROVEMENT PROGRAM 

The preceding ideas were programmed in a digital computer program 
called CIP (Compensator Improvement Program) ; A complete fortran 
version of this program is contained in the Appendix. The general 
iterating procedure employed by: CIP is as follows: 

1. Using the compensator at hand, the program calculates the 
critical points, i. e. , stability margins , attenuation margins 
and other points of interest. 

2. If this is the first iteration a preselected step size is 
chosen. Otherwise, a step size is selected according to one 
of two criteria. 

3. Next the active constraints are separated from the inactive 
constraints. 

4. After this, unit vectors in the direction of the gradients 
with respect to the variable compensator coefficients are 
obtained (The numerator partials are listed first). 

5. Then using a chosen dot product vector the , directional vector 
is determined (For the normalized gradient vectors calculated 
in 4, a suitable dot product vector has been found to be a 
vector whose components are l T s). 

•k 

The other points of interest are frequency response points on 
which special attention is to be placed, for example, points to be 
properly phased, certain gain or phase margins, etc. 
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6. Finally, the directional vector is normalized with respect 
to its magnitude; the compensator coefficients are changed 
according to the normalized directional vector and the step 
size; then, the complete process is repeated. 

In order to initiate the program, an input of discrete open loop 
frequency responses in the form of frequency and real and imaginary 
parts are required. Allowances are made for five channels of such 
information with a maximum of 999 points for each channel. This means 
that in Step 1 the actual critical points of the frequency response 
are not located — only approximate values are found. However, exper- 
ience has shown that the approximate values suffice. 

In order to determine better approximations to the critical 
points the input would require open loop transfer functions (Equation 
5-3) for each channel. The more accurate approximations of the 
critical points could be found by finding the real roots of equations 
of degree 2n, where n is the total number of the open loop system (See 
5-1).* For systems above tenth order this is completely impractical 
due to the amount of computation time necessary to perform this task. 
Furthermore, in many practical situations an experimental discrete 
frequency response is the best information available for describing 
the system. In other words an experimental frequency response is 
obtained, and using this data a transfer function of the system is 
approximated . 

Also, some initial compensator for each allowable channel is 
required. The amount of initial compensation must be enough to 

*In this discussion it is assummed that due to round-off error 
a computer is not capable of ' getting exact; solutions ‘ of non- integer 
problems. 
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stabilize the system.* If the system is open loop stable then each 
initial compensator can be chosen as an equivalent 1 compensator, i.e., 
the numerator and denominator factors are chosen to be. the same. The 
compensators may be either in a factored or unfactored form (It is 
apparent that the unfactored form is just a special case of the 
factored form) . 

In Step 2 the proper step size is chosen. In the CIP one. of two 
procedures for selecting the step size is employed. These are 

a. Require the betterment of all active constraints from the 
last iteration. 

b. Require the sum of the differences of all active constraint 
values and their desired values to increase from the last 
iteration (For this sum all active constraints of the <_ 
form have been changed- to the >_ form by multiplying by -1). 

Procedure a indicates the program is to be used in the TIFR phase, 
while procedure b designates the program as SIFR. The choice of the 
criteria used is left to the designer. If the one chosen is satisfied, 
the present step size is doubled, provided that the doubling process 
does not exceed some . preselected maximum step size value.** Otherwise, 
the maximum step size value is utilized. Regardless of which of these 
occurs the program continues to the next iteration. On the other 
hand, if the continuance criterion is not satisfied then the step size 
is halved and the present iteration is repeated if the step size is 

k 

If the system is not stable then relative stability has no 
meaning — although relative instability might. 

k k 

The main reason for limiting the step size is to keep the com- 
pensator from becoming unstable on a single iteration. 
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greater than some chosen minimum step size. When the step size 

* 

becomes less than the minimum, value the program is terminated. 

Steps 3, 4, and 5 are simply operations necessary for employ- 
ment of the constraint improvement algorithm of Chapter V, whereas, 
in Step 6, the compensator coefficients are actually changed. In 
Step 5 the reason for reducing the directional vector to a unit 
vector is so that the step size actually designates the overall change 
in the compensator coefficients. Otherwise this would not be the 
case, 

The output. of the CIP can be controlled to occur at every 
iteration or at set increments, i. e. , a set number of iterations 
can be skipped between outputs. At any iteration at which an output 
occurs the following information is printed by the CIP: 

1. Iteration number 

2. Constraint values 

3. Frequencies where constraints occur 

4. Desired constraint values 

5. Type of . constraints 

6. Directional vector at the last iteration 

7. Compensators at the present iteration 

In 5 the type of constraints denotes whether it is a phase margin, a 
gain margin, a stability margin, or an attenuation margin, and the 
symbols used to denote these are respectively P, G, S, and A. 

In the program stability margins are the main vehicles for 
determining the relative stability of the system. The concepts of 

* 

The program, also, has a maximum iteration termination condition. 
Since this has no effect on convergence, it was not included. 
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classical phase and gain margins have been included in the program 
because in some special cases these can be used to control proper 
phasing and various dynamical responses of the system. Furthermore, 
it should be pointed out that the measurement of these concepts is 
carried out exactly as stability margins, i. e. , distances from the 
-1 + jO point. Of course there is a one-to-one correspondence 
between this measuring method and the normal methods of measuring 
phase and gain margins. 



CHAPTER VII 


LARGE SYSTEM EXAMPLES 

In order to illustrate the practical usefulness of CIP, the 
improvements . of the compensators for large systems are presented. 

This is done by way of two examples; the first example is a single 
channel system, while the second example is a dual channel system. 

The two systems are not the same , although they are very similar. 

Single Channel Example 

In this example the system under consideration is similar to 
that shown in Figure 1, but only one channel is fed back. The 
system’s dynamics, 0i(s)/R(s), are described by the gain vs frequency 
and the phase vs frequency plots shown in Figures 4 and 5. This sys- 
tem is a model of the Saturn V/Sl-C Dry Work Shop at a flight time of 
80 seconds. By an inspection of these frequency response plots it 
is revealed that this system has several poles near the jw-axis. 

This deduction is based on the spike shaped gain response and the 

i 

almost discontinuous changes in the phase response. These poles near 
the j w-axis are due to various sloshing and bending modes of the 
vehicle. 

This vehicle is inherently open loop unstable. Thus, it is 
necessary to use a control scheme, such as depicted by Figure 1, to 
stabilize it. Also, unity feedback with a pure gain compensator is 
not sufficient to stabilize the system. A compensator with unity 
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Frequency in CPS 

Figure 4. Gain vs Frequency for Uncompensated System 
of Single Channel Example 
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Frequency in CPS 


Figure 5. Phase vs Frequency for Uncompensated System 
of Single Channel Example 
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feedback which is capable of stabilizing the system is 


1.0 + 11.79440s + 28.59200s 2 lOOvO + 6.05720s + 7.56640s 2 

1.0 + 21.56500s +6. 05650s 2 100.0 + 10.06500s + 6.32880s 2 


1000.0 +. 19.08700s + 3.73500s 2 

1000.0 +330. 35200s + 19.02000s 2 


(El-1) 


The GH(joj) compensated frequency response is shown in Figure 6. In- 
cluding the compensator, this frequency response represents a 29th 
over a 35th order system. 

In the design of the preceding compensator several physical 
limitations and constraints were considered — other than just stability 
of the system (In fact, stabilization of the system can be easily 
accomplished by a simple lead network with a reduced d. c. gain). 

Some of these are 

1. From past history it is known that compensators with very 
small d. c. gains produce poor wind responses. An acceptable 
value of d. c. gain is 0.9. 

2. On the GH(ja)) frequency response the first negative real axis 
crossing with respect to increasing frequency is called the 
aerodynamical gain margin. Experimentation has shown that 
the major effect of an "engine^-out" is a reduction of this 
margin. A safe crossing point is considered as -2 or less 
(or a frequency response magnitude greater than 2). 

3. For a small band of frequencies around 1.199 Hz the frequency 
response is dominated by the first bending mode. It is 
desirable to attenuate this band of frequencies. However, 

to even approach other system requirements and perform this 
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90 ° 



270° 


Figure 6. Initial Compensated GH(jw) Frequency Response 
for the Single Channel Example 
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attenuation has been* practically impossible. It has been 
found that the same effect results if this band of frequen- 
cies is phased in the’ right half of the GH(joj) plane. Due 
to the fact that the frequency of ; this mode is not known 
exactly, it is necessary- to require- larger phase margins for 
this mode than normally required. Acceptable margins are a 
lead phase margin of about 55° and a lag phase margin of 
about 90° (The reason for the difference is that in most 
physical systems phase lag is more probable to occur than 
phase lead) . 

4. For frequencies greater than 2.1 Hz the GH(jto) frequency 
response is dominated by the higher order bending modes. The 
control system can be deterred from resonating at any of 
these higher modes by attenuating to a certain degree all 
frequencies above 2.1 Hz. These frequencies are considered 
satisfactorily attenuated' if the magnitude of the GH(jw) 
frequency response is less than 0.25 for f > 2.1 Hz. 

5. Besides the above frequency response requirements, it is 
desirable for all' stability margins to be 0. 5 or greater 
(Notice that in terms of classical stability margins this is 
approximately equivalent to having phase margins of 30° and 
gain margins of 2 or better). 

By an observation of Figure 6 it becomes evident that all of the above 
specifications are not met. This becomes even more obvious after an 
inspection of Table 1. In this table the first margin is the aero- 
dynamical gain margin and the next two margins are the lead and lag 
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k 

phase margins of the 1st bending mode, respectively. The remaining 
margins listed under attenuated frequency information are stability 
margins as defined in this paper, and of course the attenuated infor- 
mation is representative of the attenuation margins above f = 2.1 Hz. 

In the CIP program the following specif ications were made: 

1. Determine the aerodynamical gain margin and improve it if 
it is less than 2. In order to improve any point it is 
necessary to specify what point or points in the complex 
GH(ja)) plane this point'is to be pulled or pushed with 
respect to. For this example it is chosen to push this 
point with respect to the -1 + jO point. 

2. Determine the lead and lag phase margins of the first bend- 
ing mode and improve either or both if they fall below 0.9 
and 1.3, respectively. To improve these it is chosen to 
push them from the -1 + jO point. 

3. Detect all stability margins and increase those less than 
0.505. Again the -1 + jO point is chosen as a pushing point. 

4. Detect all attenuation margins for f > 2.1 Hz. and decrease 

all of those greater than 0.25. For these margins the origin of 
the GH(jo)) plane is chosen as a pulling point. 


k 

The measurements of these stability margins are made in the 
same manner as stability margins defined in Chapter II, i.e. , the 
distance from the -1+jO point. Measuring gain margins in this way 
is quite natural. However, measuring phase margins ' in this way is not 
as straight forward, even though there is a one to one correspondence. 
The equations relating the two are: d = 2 sin 6/2 and 0 = 2 arcsin d/2, 
where d is the distance from the -1 + jO point and 0 is the phase 
margin. Of course d is limited to the closed interval [0,2] . 
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The continuance criterion chosen was b of Chapter VI. With these 
insertions and the necessary frequency response information in CIP, 
the following compensator was obtained after 2000 iterations (or 
approximately 30 minutes on a UNIVAC 1106) : 


1.0 + 74.40524s + 107.13383s 2 100.0 + 7.29719s + 8.68710s 2 
c Cs ' * y 1.0 + 124.68711s + 16.85849s 2 100.0 + 11.98668s + 9.15484s 2 


1000.0 + 12.10541s + 3.11162s 2 

1000.0 + 219.54201s + 20.42297s 2 


(El-2) 


A tableau of the pertinent information at iteration 2000 is shown in 
Table 2. From this tableau it is seen that most margins are, for 
practical purposes, satisfied. The reason that several of the margins 
have values that are only approximately equal to the desired values is 
that, in most instances, after a margin becomes inactive it has a 
tendency to oscillate between activity and inactivity on higher itera- 
tions, However, by establishing an upper limit on the step size from 
iteration to iteration these constraints are coerced to remain in a 
vicinity of their desired values (For this example the maximum step 
size was chosen as 0,1 for the first 1000 iterations; then, to speed 
up convergence it was changed to O'. 2 for the next 1000 iterations). 

The three smallest stability margins do not belong in the same 
category as those mentioned above "because' at' no time were they inactive. 
Since program termination was maximum iterations, additional improve- 
ments in these constraints is quite conceivable. Nevertheless, the 
convergence curve shown in Figure 7 indicates many more iterations will 
be required before any appreciable change in the smallest of these 
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Figure 7. Convergence Curve for Single Channel Example 
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margins is recorded. With an occurrence such as this the designer is 
left with three alternatives: 

1. Accept the present design. 

2. Pay the toll of additional computer time and attempt 
additional iterations. 

3. Change some of the desired constraints and continue the 
program. 

From experience it has been found" that small changes in the desired 

: 

margins can result in marked effects; r ‘ : As for ~ the case 'under dis- 
cussion the GH(jo)) frequency response in Figure 8 reveals that for 
practical purposes the compensator for 'iteration 2000 ' is satisfactory . 1 0 

Dual Channel Example 

Again reference is made to Figure 1, except in this case it is 
assumed that j = 2, i. e two ' channels are fed back. The uncompen- 
sated open loop system is described by the gain and phase frequency 
responses shown in Figures 9 , 10, 11, and 12. Figures 9 and 10 
represent the gain and phaseplots of 0^ (s) /R(s) , while Figures 11 
and 12 are the gain and phase plots . of 0£ (s) /R(s) . This system is 
typical of the Saturn V/Sl-C Sky Lab at a flight time of 105 seconds . 


* 

It should be noted that at the end of iteration 2000 the C IP was 
slightly modified so that a better calculation of the first negative 
real axis crossing frequency was obtained. After this, additional 
iterations were attempted and in less- than 50 iterations ' the smallest 
stability margin was increased fromO; 46513 to 0.48177. In another 
instance the compensator whose- smallest" stability margin was 0.48177 
was used as the starting compensator in another run in which the rela- 
tive stability requirements were lowerecj to* 0.49 while the other 
system requirements were the same as’ previously stated. In less than 
50 iterations all system requirements were completely satisfied. 
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90 ° 



270° 


Figure 8. GH(ju)) Compensated Frequency Response at 

Iteration 2000 for the Single Channel Example 
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Figure 10* Phase vs Frequency for Channel 1 of Uncom- 
pensated System of Dual Channel Example 
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Figure 11. Gain vs Frequency for Channel 2 of Uncom- 
pensated System of Dual Channel Example 
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F 5 C'S UF NT Y IN M7 


Figure 12. Phase vs Frequency for Channel 2 of Uncom- 
pensated System of Dual Channel Example 
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Compensators which have. been r designed for this system are 


„ , \ . 1000.0 + 6.54732s +^4 • 

G. (s) = 1.26 , 

1 100.0 + 1.43424s 


100.0 + 6.04029s 

100.0 + 6.17455s 


10.0 + 3.69000s 0.1 +, 1.04000s . Iv0_ 

10.0 + 2.32980s 0.1 + 2.33536s 10.0 + 1.05603s 

100.0 

100.0 + 4.13275s 


( s ) 


1000.0 + 2.91040s + 4.50787s 2 100.0 + 4.71096s 


100.0 + 3.52502s 


100.0 + 4.61899s 


10.0 10.0 

100.0 + 5.49396s 10.0 + 1.21426s 


10.0 + 2.85080s 


With these compensators inserted in the system the compensated open 
loop GH(joj) frequency response, C (j to) /R(ju)) , with the loop broken at 
a is that shown in Figure 13. 

It is desired to make several improvements in this frequency 
response. These conditional improvements are 

1. Keep the aerodynamical gain margin at 4.37 or greater. 

2. Increase all stability margins of 0.49 or less. 

3. Maintain the lead and lag phase margins of the first 
bending mode at 55° and 90° or better. 

4. Decrease all attenuation margins occurring at frequencies 
above 2.0 hz when 0.2 or greater. 


^c. 
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In order to make these improvements- the following specifications are 
made in CIPs 

1. Whenever the aerodynamical gain' margin is 4.8 or less it is 
pulled with respect to-’ the' -7 -j 3 point and pushed with 
respect. to the -1 +j0 point. 

2. All stability margin points less than 0.49 are pushed with 
respect to the -1 + jO. 

3. The lead and lag-phase margins are pulled with' respect to 
the 1 + jO point when less than 0; 9 and 1.3 respectively. 
Also, the attenuation margins' occurring at" frequencies 
between these two are decreased by pulling with respect to 
the origin of the GH(j to) plane if ; they are greater than 9.0. 

4. The attenuation margins above 2.0 hz are : decreased by pulling 
them with respect to the origin. 

With these specifications, 357' frequency response points for each 
channel, and the initial compensators, ^E2-l) and (E2-2), in the C IP, 
the following compensators were obtained after 200 iterations or about 
10 minutes on a Univac,1106: 

_ , . , 1000.0 + • 7i 07293s + 7 v 02583s 2 • 100.0: 

G (s) = 1.26 

2 100.0 + 1.21230s 100.0 + 10.48567s 

10.0 + 3.43938s - Ovl + -1,21370s 

10.0 + 1.14372s 0.1 + 2.51497s 

1.0 . 100 . 0 - 


10.0 + 1.14372s 100.0 + 9.34985s 


(E2-3) 
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1000.0 + 6.74527s + 4.53868s 2 100.0 + 4.57638s 

100.0 + 0.0s 100.0 + 1.26840s 


10.0 10.0 

100.0 + 6.96980s= 10.0 + 1.44585s 


10.0 

10.0 + 1.44585s 


(E2-4) 


An evaluation of the amount of improvement can be made by comparing 
the initial tableau, Table 3, of important information to the final 
tableau, Table 4. As in the last example the first margin is the 
aerodynamical gain margin, and the next two margins are the lead and 
lag phase margins of the first bending mode respectively. The remain- 
ing margins under relative stability information are listed as stability 
margins. The margins under the attenuated frequency information are 
the attenuation margins above 1.2 hz. The desired margins' values are 
listed in the right hand column. 

Taking into account the desired improvements it is seen that 

t 

significant improvement has been made. Furthermore, this is reinforced 
by comparing the initial compensated frequency response, Figure 13, to 
the compensated frequency response at iteration 200, Figure 14. The 
termination reason was maximum iterations; thus, as in the first 
example the designer is left with the same three alternatives. From 
the convergence curve shown in Figure 15 , it appears that several 
additional iterations may have to be attempted before any significant 
improvement in the smallest stability margin is observed. The impor- 
tance of this example is the significant improvement over the initial 
frequency response. 
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Table 3. Initial Table of Values for Dual Channel Example 
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Figure 15. Convergence Curve of Dual Channel Example 
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Additional Analysis of Results and Comments 

The results obtained from the two examples clearly indicate that 
the CIP can be a valuable design aid. It must be pointed out that as 
the name, Compensator Improvement 'Program, imples the program is a 
design aid, not a design technique; ' That is, the program does not 
decide the order, the type, or the number of compensators necessary. 
All of this requires good engineering- judgement before the running of 
the program is attempted. 

As the two examples exemplified - the ' solution cannot be worse than 
the original compensator if the specifications- on the input are made 
properly. In regard to stability margins and attenuation margins this 
simply requires pushing and pulling these, respectively, with respect 
to the -1 + jO and 0 + jo points. By doing this, these can always be 
bettered, except when they proceed from activity to inactivity. How- 
ever, the amount of slippage in going from inactivity to activity can 
be minimized by choosing a reasonable maximum step size such as 0.1 or 
less of the smallest compensator -coefficient. As long as a margin 
stays in a vicinity of the desired value it is acceptable. 

The specif ications for insuring the ' improvement in gain and 
phase margins are not always' as"simple ' as those for" stability margins 
and attenuation margins. In fact; in : many instances 'it is necessary 
to push and pull these with respect" to • two points in the" complex 
plane. This is especially true ' if ; the acute 'angle between the tangent 
to the GH(jto) frequency response where ' these occur and either the tan- 
gent to the unit circle or real azis is very small. Both of these 
cases are illustrated in Figure 16 where tangents to some hypothetical 
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GH(joj) frequency response are assumed as A and B. The points a and $ 
are the points where the margins occur. If they are perturbed so that 
the distances between them ■ and the -1 +• jO point are increased, then 
they are allowed to move in any direction Which has a positive dot 
product with vectors emanating from the -^1 +:j0 point to these points. 
Suppose that a was perturbed ' in the direction 0 ' indicated ' in Figure 16. 
It is obvious that by moving' a in this - direction the vector from -1 + 
jO to a is increasing in magnitude; ' However ,' after a is perturbed it 
is no longer the point of interest; Some other point such as A is 
then the point under consideration, where : A is- in some neighborhood of 
a. From practical considerations it' is known that if a moves in the 
direction 0 then a small neighborhood around a will move : in the direc- 
tion 0. Let A be in this neighborhood; "The result is that A will be 
the new point of intersection with the real axis, and, furthermore, 
its distance from the -1 + jO point is less than what a's was. Similar 
results can be demonstrated for B. 

These types of problems can be' circumvented by perturbing a point 
with respect to two points in the complex plane. In fact consider the 
example in the last paragraph. "' Suppose that a is not only pushed with 
respect to the -1 + jO point, but it is also pulled with respect to 
the -7 - j4 point. The permissible region for the movement' of ' a now 
becomes the intersection of the permissible region for pushing from 
the -1 + jO point and the permissible region for-pulling with respect 
to the -7 - j4 point. The result is the cross-hatched area in 
Figure 16. Movement of a anywhere in this region cannot result in 
the gain margin being decreased. 
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For the single channel example" conditions did not exist to warrant 
pertubations with respect to more than one point.. On the other hand 
the dual channel example required 'perturbing the aerodynamical gain 
margin with respect to two points; ' Runs ' in which this was not done 
resulted in a significant reduction in this margin. 

In neither example did the' lead : and lag -phase margins of the first 
bending mode become active. In the single" channel example; conditions 
just never prevailed. As for ■ the- dual channel example, conditions . 
would have probably resultedif the magnitude of r the 'first bending mode 
had not been controlled by the attenuation-margin technique; Since the 
frequencies where these margins ' occur" are' very close ' to the frequency 
of the first bending mode, then it ’ is quite : natural that an increase : 
in the first bending mode magnitude would have resulted in the reduc- 
tion of at least one of these margins. 

The program indicatedfor the dual' channel example that better 
results could be obtained with one less zeroin the numerator of the 
first channel's compensator and one less pole in the second channel's 
compensator. It did this by driving these" to infinity. It also drove 
two poles in each channel to equal values. This probably indicates 
that if these poles were included in second order factors they would 
split into complex conjugates. ' "However , the first order pole : factors 
were chosen so that complex poles "would not be allowed. 

One other fact which should be pointed out : is that the program 
was used in the SIFR mode. However ; "because of the maximum step size 
choices (0.1 for the first 1000' iterations of the first example and 



0.2 for the second 1000 iterations and 0.1 for the second example) the 
program actually performed in the TIFR mode.* 

One phenomenon which ‘ should' not' pass'without'mention is the 
apparent unsmoothness of the convergence curves ; Figures 7 and 15. In 
actuality, these curves should be discrete' curves* For convenience 
they were drawn continuously; The "sharp',' abrupt changes, where the 
smallest stability margins make much greater gains thanonother itera- 
tions, occur at iterations where- the aerodynamical gain margin became 
inactive. This allowed the' smallest 'stability margin to- make a marked 
gain for one iteration; While- this was' occurring ' the aerodynamical 
gain margin was returning ' to activity; " * Once" it" became" active ' again 
the rate of increase- of ; the' smallest" stability-margin decreased. On 
higher iterations the curve was smooth" until the ^aerodynamical gain 
margin went inactive" again ; ' at' which' time 'the "process was' repeated. 

The overall effect of the program is a rU ratchet" r type, i; e. , once a 
margin' is increased; it will not decrease. 


Of course again this is neglecting instances where constraints 
went from inactivity to activity. 



CHAPTER VIII 


CONCLUSION, LIMITATIONS , AND SUGGESTED 
FUTURE STUDIES 

Summary 

In this dissertation, the theory for a compensator improvement 
algorithm has been presented; • The r goal . from - the - onset was to accom- 
plish this by way of mathematical' programming. "Thus, in Chapter I 
a concise review of the more popular mathematical programming tech- 
niques was given. After this - review a discussion of the uses of 
mathematical programming in the' design of "control systems was pre-r 
sented. Also, a discussion of the uses of mathematical programming 
in the design of control systems was made, * In this discussion it was 
pointed out that only a small ■ amount ■ of * effort' has - been devoted to 
using mathematical programming- as* an* aid' in the design of control 
systems by classical means. Furthermore, 'it was' shown that the tech- 
niques which have been developed suffer- from' some" serious drawbacks . 
Thus, the thesis of this' dissertation" was' to develop a- computerized 
compensator - design* procedure" which* circumvented these pitfalls . 

In Chapter' Hi some important" concepts for the*measuring of 

expected performance' of - a - control system were* given.' • This " involved 
defining relative" stability in* a 'way* somewhat" different from the 
normal textbook definition. "Alsoi -concepts of relative attenuation 

and proper* phasing were defined; "Finally; using" these the design of 
a compensator was formulated as a mathematical programming problem — 
which in the end resulted in' a' strict" constraint problem. 
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In Chapter III compensator - limitations for two possible'iterative 
techniques for solving the problem- formulated in Chapter II were pre- 
sented by the proving of two - theorems ; ‘ The first" theorem showed that 
to be assured of being’ able'to'perturb'n'points'in* the GH(jw) plane 
in n optimal directions ' there' must ’exist' 2n coefficients ’ for variance. 
On the other hand; Theorem'2 stated thafif "each point 'was given 180° 
of freedom for movement (a sub-optimal” direction); then only n coef- 
ficients were ’needed' for'variance; ' From this 'it "was deduced that a 
sub-optimal algorithm" would be * the' most” practical . 

Then, in Chapter IV the" development' of' a sub-optimal ’algorithm 
was made. The* result was' the' evolvement of " the’ constraint improve- 
ment algorithm. In this development" several definitions were given, 
e.g., total improved frequency response, sum improved frequency 
response, improved solution, and active and passive constraints. 

In order to employ the constraint improvement algorithm in 
Chapter IV, it was expedient to' have the gradients” of the active con- 
straints. These were found in Chapter V for a general j*"* 1 channel 
control system. Furthermore, the’partials were derived so that” push- 
ing or pulling on points of the' frequency response* could be accom- 
plished with respect to any points* desired'in the "complex ”GH(j(jj) 
plane. 

Next , the ideas and material ’ in Chapters II; III; IV, and’ V were 
included in a computer program' called C IP ' (Compensator Improvement 
Program). In Chapter VI the general iterating procedure of this pro- 
gram was incorporated. In addition; several special programming 
techniques employed by CIP were presented in this chapter. 
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Chapter VII was used to demonstrate the practicality of CIP. This 
was illustrated by two large system examples. These examples clearly 
showed the program's capability of handling single or multi-channel con- 
trol systems. A significant amount' of -improvement in the frequency 
response of both systems was' seen after an application of CIP. Also, 
curves to show the convergence properties of "GIP were given. In 
addition, several comments in "regard" to proper specifications for the 
program were mentioned. 

Limitations and Concluding Remarks 

One of the limitations of CIP is' that the initial compensator must 
be chosen to stabilize the system." This" is' the reason' that the program 
was termed an "improvement program 1 ’ rather than a design program. A 
major goal of the program is to improve stability margins, etc., from 
one iteration to another. ' Obviously ; ' if - the system is initially un- 
stable then stability margins' have no meaning. 

Another shortcoming of CIP" is that a choice' of the components of 
the c vector in Chapter IV must'be'made. ' If the strict constraint 
problem has a solution which is reachable from the initial starting 
point, the choice of the c vector- has ■ little consequence other than to 
affect the rate of convergence; " 'However;' "if the problem does not have 
an obtainable solution, then the choice of this vector will definitely 
determine the relative extremal' where - convergence occurs : ' Neverthe- 
less, it should be pointed out that - if the initial guess at the 
solution is not a relative extremal then the - solution' at convergence 


will be better than the initial solution. 
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A very good property which CIP possesses is an inherent ability 
not to design an unstable' compensator; provided the step size is main- 
tained at a reasonable value; The’ reason for this is that CIP con- 
tinuously improves relative - stability; • thus - the stability of the 
system cannot decrease. 

Although CIP requires 'a' choice' of - the". c vector- elements, it still 
has the capabilities of yielding'a'practical'design'on every run. As 
long as the input specifications" of the program- are" properly made, 
the program cannot yield a compensator worse than the original compen- 
sator. CIP is not a design' technique; but it is a design aid. 

Suggested Future Studies 

There are several areas in which the work in this dissertation 
can be extended. One such study could involve using the constraint 
improvement algorithm in other design problems in engineering and 
science. This author does not see -any reason that it could not be 
used to make improvements in any design -where the number of -variables 
is greater than the number of constraints to be controlled and where 
the gradient vectors of the constraints are deterministic. 

Also, it is foreseen by this author - that the -constraint improve- 
ment algorithm could be the basis of -a new or extended gradient algo- 
rithm for nonlinear programming; 'For example, if any - of the elements 
of the c vector are set to zero then 'the "determined directional vector 
will lie in thetangent planes "of the constraints corresponding to the 
c's with zero value. Of course this would be similar to the gradient 
projection technique mentioned ‘in Chapter- I; - However; it - is deemed by 
this author that by using the constraint improvement approach an 
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optimal gradient projection algorithm can be developed. Up to the 
present such an algorithm has not been developed. 

In regard to future studies in compensator design, there would be 
nothing wrong with starting with physical electrical networks, rather 
than transfer functions. If a program started with a network and varied 
the elements for making the improvements described previously, the end 
results would be the actual network needed. The practicality of this 
network would depend upon the constraints placed on the network ele- 
ments. 

A compensator design procedure could be devised using the con- 
straint improvement algorithm on the Routh-Hurwitz array. By forming 
the characteristic equation as a function of the compensator coeffic- 
ients, the first two rows of the Routh array can be formulated as 
functions of these compensator coefficients. Since it is known how the 
other rows of the array are formed from the first two rows, the changes 
in the elements of the first column of the array with respect to the 
compensator coefficients could be determined by an application of the 
chain rule for partial derivatives. Then, the constraint improvement 
algorithm could be used to drive all the negative elements of the first 
column positive, as long as the number of negative elements did not 
exceed the number of compensator coefficients. If all the elements are 
driven positive then a certain amount of relative stability could be 
achieved by evaluating the characteristic equation at (s + a) where a is 
a positive real number; the previously mentioned procedure can now be 
applied to the new characteristic equation. If in this application all 
elements of the first column could again be driven to positive values, 
then it would be known that no pole of the closed loop system has a 
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real part greater than - a. This process could be repeated until a 
desired value of a is achieved or until all the elements of the first 
column of one of the characteristic equations cannot be driven 
positive. 



APPENDIX 


COMPENSATOR IMPROVEMENT PROGRAM 

In the following is a complete Fortran version of the Compensator 
Improvement Program. The program is completely self-contained, i.e., 
it does not require any system library, etc. The necessary input to 
the program is explained in the comment statements at the beginning of 
the main program. Furthermore, all inputs except the frequency 
response points are printed out with explanations of the input speci- 
fications. The other output is, also, explained by certain comments 
printed out with the information. 



oooooonoononoooonnoonnnooononoooooonoooonoonoonooo 


83 


MAIN PROGRAM 

DEFINITIONS OF I/O VARIABLES 

KCHNL -NO. OF CHANNELS FED BACK 

NUMC(I) -NO. OF COMPENSATORS IN I-TH CHANNEL 

NRATOR ( I » J) -NUMERATOR ORDER OF J-TH COMPENSATOR IN THE I-TH CHANNEL 
NDENOM( I » J) -DENOMINATOR ORDER OF J-TH COMP. IN I-TH CHANNEL 
XCOMN(I.J) -NUMERATOR COEFFICIENTS OF J-TH COMP. IN I-TH CHNL. 
YCOMN( I » J) -DENOM. COEFFICIENTS OF J-TH COMP. IN I-TH CHNL. 

OMEGA ( I ) -I-TH FREQ. (ASSUMED TO BE IN HZ.) 

GRA ( I » J) -J-TH REAL PART OF OPEN LOOP FREQ. RESP. OF I-TH CHNL. 

GIA(I»J) -J-TH IMAG. PART OF OPEN LOOP FREQ. RESP. OF I-TH CHNL. 

KSTART -STARTING ITERATION NO. 

KQUIT -STOPPING ITERATION NO. 

KPOINT -NO. OF POINTS FROM OPEN LOOP FREQ. RESPONSE USED 
KPRINT - NO. OF ITERATIONS SKIPPED BETWEEN PRINTING OF INFOR. 
STPMAX -MAXIMUM CHANGE TO BE MADE IN COMPENSATOR COEFFICIENTS 
ON ANY ONE ITERATION (PROBABLY NO MORE THAN 30* OF THE 
SMALLEST COMPENSATOR COEFFICIENT OF THE INITIAL 
COMPENSATOR) 

STPMIN - MINIMUM STEP SIZE DESIGNATOR 
F10 & FH - FREQUENCIES BETWEEN WHICH G.M.'S ARE FOUND 

F12 & F13 - FREQUENCIES BETWEEN WHICH P.M.'S ARE FOUND 

FMIN - A.M.'S ARE FOUND FOR FREQS. ABOVE THIS FREQ. 

variables for gain margin radii designations 

IF FREQ. .LE. FI DESIRED MARGIN = R1 

IF FREQ. .Gt. FI BUT .LT. F2 DESIRED MARGIN = R2 

IF FREQ. .GE. F2 DESIRED MARGIN = R3 

variables for phase margin radii designations 

IF FREQ. .LE. F3 DESIRED MARGIN = R4 

IF FREQ. .GT. F3 BUT .LT. F4 DESIRED MARGIN = R5 

IF FREQ. .GE. F4 DESIRED MARGIN = R6 

VARIABLES For STABILITY margin radii designations 
IF FREQ. .LE. F5 DESIREO MARGIN = R7 

IF FREQ. .GT. F5 BUT .LT. F6 DESIRED MARGIN = R8 

IF FREQ. .GE. F6 DESIRED MARGIN = R9 

VARIABLES For attenuation margin radii designations 
IF FREQ. .LE. F7 desired MARGIN = RIO 

IF FREQ. .GT. F7 BUT .LT. F8 DESIREO MARGIN = Rll 

IF FREQ. .GE. F8 DESIRED MARGIN = R12 

GAIN(I)-DEN0TES INITIAL D. C. GAIN VALUE FOR I-TH CHANNEL 

KNR(I) -NUMBER OF NUMERATOR COEFS. FOR 1-TH CHANNEL 
KDR(I) -NUMBER OF DENOM. COEFS. IN I-TH CHANNEL 
KONT ( I )-D.C . DESIGNATOR FOR I-TH CHANNEL 

KONT(I)=l GAIN ALLOWEO TO VARY 
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K0NT(I>=2 GAIN NOT ALLOWED TO VARY 

KIFM -NO. CHANNELS THAT FREQ. RESP. INFIRMATION IS TO BE READ IN 

PPT(I) -POINTS THAT THE CRITICAL FREQUENCIES WILL BE 
PERTURBED WITH RESPECT TO (COMPLEX POINTS) 

1=1 GAIN MARGIN POINT 
1=2 PHASE MARGIN POINT 
1=3 STABILITY MARGIN POINT 
1=4 ATTENUATION MARGIN POINT 

lsn(D - denotes whether points are to be pushed or pulled 
lsn=-i point to be pulled 

LSN=+1 POINT TO BE PUSHED 

INCGMS -INDICATES WHETHER G.M.'S ARE TO BE ARTIFICALLY 
INCLUDED AS S.M.*S 
INCGMS=0 NOT INCLUDED 

INCGMS=1 INCLUDED 

INCPMS -INDICATES WHETHER P.M.*S ARE TO BE ARTIFICALLY 
INCLUDEED AS S,M.*S 
INCPMS=0 NO INCLUDED 

INCPMS=1 INCLUDED 

some interior variable definitions 


GCR(IrJ) 
GCI ( I * J) 
GCOMR ( I * J) 
GCOMR ( I * U) 
GR(I) 

GI(I) 


-REAL PART OF COMPENSATOR FREQ. RESP. AT SOME ITERATION 
-1MAG. PART OF COMPENSATOR FREQ. RESP. AT SOME ITERATION 
-REAL PARTS OF I-TH CHNL. OPEN LOOP FREQ. RESP, 

-IMAG. PARTS OF I-TH CHnL. OPEN LOOP FREQ. RESP. 

-REAL PARTS OF TOTAL OPEN LOOP FREQ. RESP. 

-IMAG. PARTS OF TOTAL OPEN LOOP FREQ. RESP. 


***** there are 13 read statements ***** 

DIMENSION XCOMN(10f50)»YCOMN(10»50)»PRY(50) »PRX(50) *STBM(99) * 

1 PX(50) *PY ( 50 ) * RQ ( 99) * GR ( 999 ) » G I < 999) * OMEGA ( 999) * GRA ( 5 * 999 ) » 

2 GIA (5.999) *G(20*99) »DV (50) » WEIGHT (50) *BCOMN(10*50) * 

3 6C0MD (10.50) *GCR( 5*999) *GCI (5*999) * GCOMR ( 5 * 999) * 

4 GCOMI (5*999) * NUMC ( 20 ) * NRATOR ( 10 * 20 ) . NDENOM ( 10 * 20 ) * CNUM ( 10 ) * 

5 CDOM(IO) *KNR(10) * KDR ( 1 0 ) * COTN ( 10 * 50 ) *COTD(10*50) 

DIMENSION KACT(99) *SML(99) 

DOUBLE PRECISION G*pV*WElGHT 
DIMENSION KoNT (20 ) * KPTS(99)* GAIN(IO) 

dimension Type ( 99 ) 

DIMENSION PPT (4 ) * LSN ( 4 ) 

COMMON TYPE 

INTEGER TYPE 

COMPLEX PPT 

READ(5*5) KCHNL 

READ (5* 5) (KONT ( I ) * 1 = 1 * KCHNL) 

READ (5* 5) (NUMC(I) *1=1* KCHNL) 

WRITE (6*1) KCHNL 

1 FORMAT ( * 0 ’ * 5X * ’NUMBER OF CHANNELS FEDBACK= * , 15) 

WRITE(6*3) ( KONT (I) *1 = 1* KCHNL ) 

3 FORMAT( 'O' *5X* *D.C. GAIN CONSTRAINT DESIGNATOR FOR EACH CHANNEL ( 
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1K0NT = 1 * ALLOWED TO VARY; K0NT=2* HELD CONSTANT ) */6X» 8 ( 12* 10X) ) 
WRITE (6*4) ( NUMC ( I ) * I = 1 * KCHNL ) 

4 FORMAT C0» *5X* 'COMPENSATORS PER CHANNEL' * 1015) 

DO 2 I=1*KCHNL 

KNAT=NUMC(I) 

READ(5*5) (NRATOR(I*J) *J=1*KNAT) 

WRITE (6* 6) I * (NRATOR ( I * J) *J=1*KNAT) 

6 FORMAT ( 'O' *5X* 'CHANNEL NO. '* I2*2X* 'NUMERATOR ORDERS' *2X» 1015) 
READ(5*5) (NDENOM(I*J) * J=1*KNAT) 

WRITE (6, 7) I* (NDENOM(Ir J) *J=1*KNAT) 

7 FORMATC • *5X* 'CHANNEL NO. '* I2*2X* 'DENOMINATOR ORDERS' * 1015) 

2 CONTINUE 

5 F0RMAT(16I5) 

READ(5*10) KSTART*KQUIT*KIFM*KPOlNT*KPRINT* R1*F1*R2*F2*R3* 

1 R4*F3*R5*F4*R6* R7*F5*R8*F6*R9* R10*F7*RU*F8*R12* FmIN*F10* 

2 F11*F12*F13*STPMAX*STPMIN 

10 FORMAT(5I5/5F10.5/5F10.5/5F10.5/5F10.5/8F10.5) 

WRITE (6* 11 ) KSTART *KQUIT *KIFM*KPOINT *KPRINT 

11 FORMAT ( 'O' *1X* 'START ITeR. = ' * I5*2X* 'STOP ITER. = ' * I5*2Xf 'NO. CHNL* 
1FREQ. RESP. IN=» *I5*2X* 'NO. OF FREQ. POINTS=» * I5*2X* 'PRINT INCREME 
2NT='*I5) 

WRITE (6*25) STPM AX * STPMIN 

25 FORMAT( 'O' »5X» 'MAXIMUM DESIGNATED STEP SIZE =' *F10.5/6X* 'MINIMUM D 
1ESIGNATED STEP SIZE ='*F10.5) 

WRITE (6 * 12) 

12 FORMAT('0»5x* 'DESIRED GAIN MARGIN RADII DESIGNATIONS') 

WRITE (6* 13) F1*R1* Fl*F2*R2* F2*R3 

13 format ( ' o» *sx * » if frequency ,le. '*fio.5»5x* 'desired margin is'* 

1 F10.5/6X* ' IF FREQUENCY .GT. » *Fl0.5*2X* 'BUT .LT. ' *Fl0.5*2X* 'DESIRE 
2D MARGIN IS' *F10.5/6X* * IF FREQUENCY .GT. ' *F10.5*2X* 'DESIRED MARGI 
3N 1S'»F10.5) 

WRITE (6* 17) F10 * Fll 

17 FORMATC »*5X*'GAIN MARGINS ARE DETERMINED BETWEEN THE FREQUENCIES 
1 OF* *F10.5*2X* 'AND' *F10.5) 

WRITE (6*14) 

14 FORMAT ( * 0 * * 5X* 'DESIRED PHASE MARGIN RADII DESIGNATIONS’) 
WR1TE(6*13) F3*R4* F3*F4*R3* F4*R6 

WRITE(6*18) F12 * Fl3 

18 FORMATC ' ► 5X * ' PHASE MARGINS ARE DETERMINED BETWEEN THE FREQUENCE 
IS OF' »F10.5*2X* 'AND' *F10.5) 

WRITE (6* 15) 

15 FORMAT( 'O' »5X» 'DESIRED STABILITY MARGIN RADII DESIGNATIONS') 
WRITE(6*13> F5.R7 * F5*F6*R8* F6*R9 

WRITE (6* 16) 

16 FORMAT( »0» *5X* 'DESIRED ATTENUATION MARGIN RADII DESIGNATIONS') 
WRITE (6* 13) F7 * RIO * F7*F8*R11* F8*R12 

WRITE (6*19) FMIN 

19 FORMATC ' *5X* 'ATTENUATION MARGINS ARE FOUND FOR FREQS. ABOVE' * 

1 F10.5) 

READ(5*50) ( GAIN ( I ) *1 = 1* KCHNL ) 

WRITE (6 *20) (I*GAIN(I) *1=1* KCHNL) 

20 FORMAT ( *0' *5X*2( 'CHANNEL NO. ’* 13* IX* * INITIAL D.C. GAIN IS**F10.5* 
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1 5X) ) 

READ(5>50) ( PPT ( I ) 1 1 = 1 r 4) 

WRITE (6»22) (PPT(I) »I=1>4) 

22 FORMAT ( ’0* »5X» •PERTUBATION POINTS FOR GAIN* PHASE* STABILITY » AND 
1ATTENUATION MARGINS* RESPECTIVELY 5 * /6X* 4 ( ’REAL* *F6*2*2X* ’ IMAG. * * 

2 F6.2*2X)) 

READ (5*5) (LSN(I) *1=1*4) 

WRITE(6,23) (LSN(I) ,1=1*4) 

23 FORMAT ( *0* *5X* ’DENOTING WHETHER EACH OF THE PRECEDING POINTS ARE * 

1 * ' TO BE PUSHING OR PULLING POINTS (PUSHING=+1 * PULLING=-1)* /6X* 

2 4 ( 12 , 10X ) ) 

READ(5*5) INCGMS»INCPMS 
WR1TE(6*24) INCGMS# INCPMS 

24 FORMAT* ’O' *5X* ’DENOTING WHETHER GAIN OR PHASE MARGINS ARE ARTIFICA 

ILLY INCLUDED AS STABILITY MARGINS (NOT INCLUDED=0* INCLUDED=1 ) ’/6X, 
2 2 ( 12* 10X) ) 

KVARY=0 

DO 21 K=1 *KCHNL 
LAMP=NUMC(K) 

KNR(K)=0 
KDR(K)=0 
DO 21 I=1*LAMP 
KVARY=KVARY+NRATOR (K * I ) 

KVARY=KVARY+NDENOM ( K r I ) 

KNR(K)= KNR(K) + NRATOR(K*I) + 1 
21 KDR (K ) = KDR(K) + NDENOM (K * I ) + 1 
DO 29 1=1 f KCHNL 

29 IFIKONT(I) .EQ.1)KVARY=KVARY+1 
DO 42 K=1*KCHNL 
LNC= KNR(K) 

LDC= KDR (K) 

READ (5* 50) ( XCOMN ( K * I ) * 1 = 1 , LNC ) 

42 READ (5* 50) (YCOMN(K,I) *I=1*LDC) 

50 FORMAT ( 8F10 . 5 ) 

60 FORMAT* ’0’ *6X* ’ INITIAL COMPENSATOR COEFFICIENTS’) 

WRITE (6 r 60) 

DO 72 K=1*KCHNL 
LNC=KNR(K) 

LDC= KDR(K) 

WR1TE(6.62)K 

62 FORMAT ( ’0’ #5X» ’CHANNEL NO. ’ , I2*2X* ’COMPENSATORS - FACTORED FORM’) 
WR1TE(6*68) 

68 FORMAT* *0’ *5X* ’NUMERATOR COEFFICIENTS’) 

WRiTE(6*70> ( XCOMN ( K * I ) , 1 = 1 * LNC ) 

WRITE(6»69) 

69 FORMAT* *0’ *5X* ’DENOMINATORS COEFFICIENTS’) 

72 WRITE (6 * 70 ) ( YCOMN (K * I > * 1 = 1 * LDC ) 

70 FORMAT (» * *5X* 10F10.5) 

C MODIFICATION OF FREQ. RESP. INFOR. BY CONTANT COMPENSATOR 
DO 135 U=1*KIFM 

135 READ(5*140) (OMEGA ( I ) » GRA ( Jr I)*GIA(J*I)* 1=1 *KP0INT ) 

140 format ( 9 F 8 . 5 ) 
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IFIKIFM.GE.KCHNDGO TO 150 

K= KIFM + 1 

DO 148 J=K»KCHNL 

DO 148 I=l»KPOINT 

XV= OMEGA(I) * 6.2831853 

GRA(J»I)= -OMEGA(I) * XV * GIA(J-lrl) 

148 GlA(JrI)= OMEGA(I) * XV * GRA(J-1»I) 

150 CONTINUE 

DO 149 J=1 » KCHNL 

DO 149 I = 1 f KPOINT 

GRA(UrI)= GRA (U» I ) * GAINIJ) 

149 GIA(U>I)= GIA(Jrl) * GAlN(J) 

190 CONTINUE 

DATA STEP » KHOP t SML2 » PSQL » S8C2/1 . 0E-02 1 0 > 0 . 0 1 1 . 0E+20 >0.0/ 
11=0 
12=0 

DO 195 K=1>KCHNL 
11= KNR(K) + II 
195 12= KDR(K) + 12 
LOX= KSTART 
200 CONTINUE 

lpresv=kvary 

NM=Q 

C EVALUATION OF VARIABLE COMPENSATOR AT CHOSEN FREQS. 

DO 210 K=1 > KPOINT 
GR (K ) =0 • 0 
GI(K)=0.0 

XV=OMEGA(K) *6.2831853 
DO 209 1=1 >KCHNL 
KCQMP= NUMC(I) 

LN0T=0 

KNOT=0 

GCR l I >K ) = 1.0 
GCl ( I »K)= 0.0 
DO 208 J=1 > KCOMP 
NTR= NRATOR(I» J)+l 
NTD= NDENOM ( I > J) +1 
DO 204 M=1>NTR 

204 CNUM(M)= XCOMN ( I > M+KNOT ) 

DO 205 M=1>NTD 

205 CDOM(M)= YCOMN(I»M+LNOT) 

knot= knot + ntr 

LNOT= LnOT + NTD 
K2= NTR-1 
K3= NTD-1 

CALL P0LFV(CNUM.K2»XV»CNR»CNI) 

CALL P0LFV(CD0M.K3»XV»CDR»CDI) 

CD= CDR**2 + CDI**2 
ACR=GCR(I»K) 

ACI= GCI(IrK) 

ACOMR= (CNR * CDR + CNI * CDI)/CD 
AC0MI=(-CNR * CDI + CNI * CDR)/CD 





GCR(I»K)= ACR * ACOMR - ACI * ACOMl 

208 GCI ( I »K)= ACR * ACOMl + ACI * ACOMR 
GCOMR(I»K)= GRA(I»K)*GCr(I»K) - GIA ( I »K ) *GCI ( I » K > 
GCOMKI»K)= GRA(I»K)*GCI(I»K) + GIA( I »K) *GCR( I >K) 
GR(K)= GR(K) + GCOMR ( I » K ) 

209 GI(K)= GI(K) + GCOMI ( I » K ) 

210 CONTINUE 

C DETERMINATION OF GAIN MARGINS POINTS BETWEEN FI AND F2 

CALL GAINMG ( GR * G I » KPOINT » NM » FI 0 t FI 1 1 KPTS » STBM t OMEGA ) 

ngms=nm 

C SETTING DESIRED STABILITY RADII OF G.M.’S 
KPm=NM+1 

IF ( NM »EQ • 0 ) GO TO 213 
DO 212 I=1»NM 
TYPE(I)= »G» 

KWHICH=KPTS(I) 

FREHZ=OmEGA(KWHICH) 

IF (FREHZ.LE.fi )RQ( I )=R1 
IF(FREHZ.GT.F1)RQ(I)=R2 
IF (FREHZ .GE.F2 ) RQ ( I ) =R3 

212 CONTINUE 

213 CONTINUE 

C DETERMINATION OF P.M. BETWEEN F3 AND F4 

CALL PHASEM ( GR » G I » KPOINT » NM t F 12 1 F 1 3 » KPTS * STbM r OMEGA ) 
IF (NM.LT »KPM)GO TO 215 
C SETTING DESIRED STABILITY RADII OF P.M.’S 
DO 214 I=KPM»NM 
TYPE(I)= »p» 

KWHICH=KPTS(I) 

FREHZ=OMEGA (KWHICH) 

IF (FREHZ .LE.F3) RQ ( I ) =R4 
IF ( FREHZ. GT .F3) RQ ( I )=R5 
IF ( FREHZ .GE.F4)RQ( I )=R6 

214 CONTINUE 
KPm=NM+1 

215 CONTINUE 

IF (NM .EQ . 0 ) GO TO 221 

klast=nm 

DO 220 I = 1 » KLAST 

IF ( ( I .LE.NGMS) .AND. ( INCGMS.EQ.l ) ) GO TO 2l9 
IF ( (I.GT.NGMS) .AND. ( INCpMS.EQ. 1 ) )GO TO 2l9 
GO TO 220 

219 KPm=KPM+1 
NM=NM+1 

KPTS(NM)=KPTS(I) 

STBM(NM)=STBM(I> 

RQ(NM)=RQ(I) 

TYPE(NM)=»S» 

220 CONTINUE 

221 CONTINUE 
KSTBM=KpM 
RPT=1.0 
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NSG=1 

FQMIN=0.0 

C DETERMINATION OF STABILITY MARGINS 

CALL SRMlNS(GR»GI»KPOINT»NM»RPT»NSG»FQMlN»KpTS»STBM»OMEGA) 

C SETTING DESIRED STABILITY MARGINS 
IF(NM.LT.KPM)GO TO 216 
DO 230 I=KPMrNM 
TYPE ( I ) = *S» 

KWHICH=KPTS(I) 

FR£HZ=OMEGA(KWHICH) 

IF(FREHZ.LE.F5)RQ(I )=R7 
IF(FREHZ.GT.F5)RQU)=R8 
IF(FREHZ.GE.F6)RQ(I)=R9 
230 CONTINUE 

C CHECKING TO SEE IF ANY P.M.*S» G.M.»S» OR S.M.’S ARE EQUAL 
C IF THERE RESULTS SOME THAT ARE EQUAL ONLY THE FIRST IS RETAINED. 

DO 228 LB=2 1 KSTBM 
DO 228 I=KSTBM»NM 

IF IKPTS (LB-1 ) .NE.KPTS ( I ) ) GO TO 228 

nm=nm-i 

DO 226 L=I*NM 
KPTS(L)= KPTS(L+1) 

STbM(L)= STBM(L+1) 

RQ(L)= RQ (L+l ) 

226 TYPE(L)= TYPE(L+1) 

228 CONTINUE 
KPM=NM+1 

216 CONTINUE 
KMIN=NM 
RPT=0.0 
NSG=-1 
FQmIN=FmIN 

C DETERMINATION OF ATTENUATION MARGINS 

CALL SRMlNS(GR»GI»KPOINT»NM*RPT»NSG»FQMIN»KpTS» STEM. OMEGA) 

C SETTING DESIRED ATTEN. MARGINS 
IF lNM.LT .KPM)GO TO 217 
DO 232 I=KPM*NM 
TYPE ( I ) = *A* 

KWHICH=KPTS(I) 

FREHZ=OMEGA(KWHICH) 

IF(FREHZ.LE.F7)RQC1)=R10 
IF(FREH2.GT.F7)RQ(I)=R11 
IF (FREHZ.GE.F8 ) RQ C I ) =Rl2 
232 CONTINUE 

217 CONTINUE 
SBC1=R1 

C DETERMINING SMALLEST STABILITY MARGINS OF PRESENT ITER. AND ALL ITER. 
SML1= 100.0 
DO 290 I=1»KMIN 
IF(STBM(I) .GT.SMLDGO To 288 
SML1= STBM(I) 

286 CONTINUE 
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IF(STBM(I) .GT.SBCDGO TO 290 
SBC1= STBM(I) 

290 CONTINUE 

IF ( SBC2 • GE • SBC 1 ) 60 TO 298 

SBC2= SBC1 

IBEST= LOX 

DO 292 K=l, KCHNL 

LNC= KNR(K) 

LDC= KDR(K) 

DO 291 I=1*LNC 

291 BCQMN(K»I)= XCOMN(K,I) 

DO 292 1=1 ,LDC 

292 BCOMD*K»I)= YCOMN*K,I> 

298 CONTINUE 

C CHECKING SATISFACTION OF SYSTEM REQUIREMENTS 
DO 320 1=1 ,NM 
PORM= 1.0 

IFII.GT.KMIN)PORM=-1.0 
310 IF ( (STBM( I ) -RQ ( I ) ) *PORM) 350 *320 >320 
320 CONTINUE 

WRITE*6,330) 

330 FORMAT (• O', 15X ******* ALL SYSTEM REQUIREMENTS HAVE BEEN MET ****** 
1) 

340 CONTINUE 

CALL OTPT1 * STBM , OMEGA , KpTS , NM , XCOMN , YCOMN , KMIN , RQ , LOX , KCHNL , NUMC * 

1 NRAT0RrNDEN0M>PRX*PRY>IlrI2) 

WRITE*6,34l) IBEST 

341 FORMAT t*0*,5X, ****** BEST COMPENSATORS WITH RESPECT TO STABILITY * 
l****«//6Xr ’OCCURRED ON STEP* , 14, 2X, ’AND THEIR COEFFICIENTS ARE!*) 

CALL MULOUT ( KCHNL , NUMC , NR ATOR » NDENOM > KNR , KDR , BCOMN * BCOMD > 
WRITE*6,345) SBC2 

345 FORMAT ( '0* *21X, 'SMALLEST STABILITY MARGIN FOR THE BEST COMPENSATOR 
/ =»»F10.8> 

WRITE*6,347) 

347 FORMAT ( *0' ,5X, ****** COMPENSATORS AND COMPENSATED FREQUENCY RESPON 
1SE AT THE LAST ITERATION PERFORMED ARE AS FOLLOWS ******) 

CALL MULOUT (KCHNL, NUMC, NR ATOR, NDENOM r KNR, KDR, XCOMN, YCOMN) 

WRITE (6, 346) 

346 FORMAT* 'O', 9X» ’COMPENSATED FREQUENCY RESPONSE'//10X, 'FREQUENCY* , 

1 2X, 'MAGNITUDE* ,3X* ’ANGLE’ ) 

DO 349 1=1 ,KPOINT 

GMTE= SQRT *GR * I ) **2 + GI*I)**2> 

AGLE= ATAN2(GI(I) ,GR*I) )*57.3 
WRITE*6,348) OMEGA(I) ,GMTE,AGLE 

348 FORMAT*' ',7X,F10.5,1X,F10.5,1X,F10.5) 

349 CONTINUE 
STOP 

350 CONTINUE 

C STEP SIZE SELECTING 

IF(L0X.EQ.KQUIT)WRITE16,351) 

351 FORMAT**0»»5X» ****** TERMINATION REASON - MAXIMUM ITERATIONS ***** 
1 *) 



91 


IF (LOX «EQ •KQUIT ) WRITE (6.400) STEP 

IFILOX.EQ.KQUITJGO TO 340 

IF(LOX.EQ.KSTART)GO TO 354 

ADD=0.0 

MAD=0 

PORM=1.0 

DO 355 1=1 *NM 

IF(I.GT,KMlN)PORM=-1.0 

IF(PORM*(STBM(I)-RQ(I) ) .GE.O.O)GO TO 355 
DO 352 J=1 * NML 

352 IFtKPTS(I) .EQ.KACT(J) )GO TO 353 

mad=mad+i 

GO TO 355 

353 CONTINUE 

C IF IT IS DESIRED TO HAVt ALL CONSTRAINTS TO BE IMPROVED AT EVERY 
C ITERATION REMOVE THE C FROM COLUMN 1 OF THE FOLLOWING CARD 
C IFlPORM*(STBM(I)-SML<I> ) .LT.-1.0E-05)GO TO 360 
ADD=ADD+PORM* (STBM ( I ) -SmL ( J ) ) 

355 CONTINUE 

IF l MAD . EQ . NML ) ADD=1 . 0 
IF(ADD.LE.O.O)GO TO 360 

354 CONTINUE 
GO TO 371 

360 STEP= STEP/ 2 • 0 

IFISTEP.LT.STPMIN )WRITE<6.365)STPMIN 

365 format (* o**5X ******* termination reason - step size is less than * 

1 *F10.5*2Xr ******* ) 

IFISTEP.LT.STPMIN )GO To 340 
LOX= LOX - 1 
GO TO 450 

371 STEP=1. 41416 * STEP 
373 CONTINUE 
SML2=SML1 

IF(STEP.GT.STPMAX)STEP= STPMAX 
C OUPUT CONTROL 

IF ( KHOP . 6T » 1 ) GO TO 410 

KH0P=KPRINT 

WRITE (6 » 400 ) STEP 

400 format ( * o * * isx» ’present step size =»*fio.7) 

CALL OTPT1 ( STBM • OMEGA , KpTS * NM » XCOMN » YCOMN , KM IN * RQ t LOX , KCHNL * NUMC * 

1 NRAT0R*NDEN0M*PRX*PRY*I1*I2) 

GO TO 420 

410 KHoP= KHOP - 1 
420 CONTINUE 

C SELECTING ACTIVE CONSTRAINTS 
K=0 

DO 411 1=1 *NM 

IF(I-1.EQ.KMIN)KMIN=K 

PORM=1.0 

IF(I.GT.KMlN)PORM=-1.0 

IF(PORM*STBM(I) .GT.PORM*RQ(I) )G0 TO 411 

K=K+1 
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KPTS(K)= KPTSU) 

TYPE (K) = TYPE ( I ) 

SML(K)= STBM(I) 

KACT(K)= KPTS(I) 

411 CONTINUE 
NM=K 

nml=nm 

C CALCULATION OF GRADIENTS OF ACTIVE CONSTRAINTS 
RPT=1.0 

CALL PARCLT(XCOMN» YCOMN,GR»GI » OMEGA »NM»NRATOR*NDENOM» 

1 KCHNL . NUMC rKONT t GCOMR » GCOM I # G * PPT * LSN » KPARC » KPTS , KnR » KDR ) 

C SET DOT PRODUCT VECTOR 
DO 422 K=1>NM 

422 WEIGHT(K>=1.0 

C CALCULUTE DIRVECTIONAL VECTOR 

lre=o 

kre=o 

423 IFINM.GT .LPRESV) WRITE (6» 415) 

IFINM.GT. LPRESV ) GO TO 340 

CALL DIRVEC (G*NM» KPARC »DV» HEIGHT) 

415 FORMATt *0* »5X» ****** TERMINATION REASON - No* OF ACTIVE CONSTRAINT 
IS IS GREATER THAN THE No. OF ALLOWABLE VARIABLES ****♦•) 

DO 426 1=1 * II 

426 PRX(I)= DVd) 

DO 427 1=1 » 12 

427 PRY (I ) = DVdl+I) 

IF (KRE.EQ.DGO TO 433 

C CKECKING POSSIBLE NEGATIVENESS OF ANY COMPENSATOR COEF. 

IF (LRE.GE. I 1+12) GO TO 433 

LRE=LRE+1 

K2=0 

K3=0 

DO 431 K=1»KCHNL 
LNC=KNR(K) 

LDC=KDR(K) 

DO 429 1=1 »LNC 
K2=K2+1 

IF(XCOMN(K»I) .GT.1.0E-05)GO TO 429 
IF(PRX(K2) .GE.O.O)GO TO 429 

lpresv=lpresv-i 

kre=i 

DO 428 J=1 » NM 

428 G(U,K2)=0.0 

429 CONTINUE 

DO 431 1=1 »LDC 
K3=K3+1 

IFlYCOMN(Kd) .GT.1.0E-OS)GO TO 431 
IF(PRY(K3) .GE. 0.0)60 TO 431 

lpresv=lpresv-i 

kre=i 

DO 430 U=1»NM 

430 e(U»Il+K3)=0.0 
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431 CONTINUE 

IFlKRE.EQ.DGO TO 423 
433 CONTINUE 
PSQ= 0.0 
DO 438 1=1 » II 
PX(I)= PRXU) 

438 PSG= PSQ + PX ( I ) **2 
DO 440 1=1 » 12 
PY ( I ) =PRY ( I ) 

440 PSQ= PSQ + PY(I)**2 
PM6= SORT (PSQ) 

PSQL= PSQ 
DEL= step/pmg 

DO 462 K=1»KCHNL 
LNC= KNR(K) 

DO 462 I = DLNC 
462 COTN(K» I )= XCOMN(K»I) 

DO 464 K=1»KCHNL 
LDC= KDR(K) 

DO 464 1=1 »LDC 

464 COTD (K » I ) = YCOMN(K#I) 

GO TO 465 

450 DEL= DEL/2.0 

DO 467 K=DKCHNL 
LNC= KNR(K) 

DO 467 I=1»LNC 

467 XCOMN (K » I ) = COTN(K»I) 

DO 468 K=1»KCHNL 
LDC= KDR(K) 

DO 468 1=1 »LDC 

468 YCOMN (K » I ) = COTD(K.I) 

465 CONTINUE 
KKK=0 

DO 470 K=1»KCHNL 
LNC= KNR(K) 

DO 470 I=1»LNC 
KKK= KKK+1 

XCOMN (K . I ) = XCOMN(K.I) + DEL * PX(KKK) 
470 IF(XCOMN(K»I) .LT. 0 . 0 ) XCOMN (K » I > =0,. 0 
KKK=Q 

DO 490 K=DKCHNL 
LDC= KDR(K) 

DO 490 1=1 » LDC 
KKK= KKK+1 

YCOMN(K» I )= YCOMN (K » I ) + DEL * PY(KKK) 
490 IF (YCOMN (K» I) .LT . 0 . 0 ) YCqMN (K » I > =0 . 0 
500 CONTINUE 

LOX= LOX + 1 
GO TO 200 
END 
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SUBROUTINE PARCLT ( XCOMN , YCOMN * GOR » GO I • OMEGA , NFREQ * NR ATOR f NDENOM » 

1 kchnl »numc»kont» gcomr » gcomi »g>ppt»lsn»nparc*kpts»knr»kdr) 

C 

c program for calculating the change of a frequency response WITH 
C RESPECT TO A CONPENSATOR COEFFICIENTS 
C 

C DEFINITIONS OF I/O VARIABLES 
C 

C XCOMN(I»U) -NUMERATOR COEFs. of compensator in i-th channel 
C YCOMN ( I » J ) -DENOM. COEFS. OF COMPENSATOR IN I-TH CHANNEL 
C GOR C I > -I-TH REAL PART OF OPEN LOOP FREQ. RESP. 

C GOKI) -I-TH IMAG. PART OF OPEN LOOP FREQ. RESP. 

C OMEGA ( I ) -I-TH FREQUENCY RESPONSE POINT 

C NFREQ -NUMBER OF MARGINS TO BE IMPROVED 

C NRATOR ( I » J ) -NUM. ORDER OF U-TH COMP. IN I-TH CHANNEL 

C NDENOM ( I .U) -DEN. ORDER OF U-TH COMP. IN I-TH CHANNEL 

C KCHNL -NUM. OF CHANNELS 

C NUMC(I) -NUM. OF COMPS. IN I-TH CHANNEL 

C KONT(I) -GAIN CONTROL NUM. FOR I-TH CHANNEL 

C GCOMR ( I > U) -REAL PART OF U-TH CHANNEL COMP. FREQ* RESP. AT U-TH FREQ. 
C GCOMI ( I > U ) -IMAG. PART OF U-TH CHNL. COMP. FrEq. RESP. AT U-TH FREQ. 

C GU»J) -U-TH PARTIAL OF I-TH FREQ. 

C Z -NEG. OF POINT FOR WHICH PARTIALS ARE DESIRED 

C L -NO. OF POINTS TO TREAT AS STABILITY MARGINS<THE REMAINING 

C ARE CONSIDERED AS ATTENUATION MARGINS) 

DIMENSION C(10) rD(10) »E(10) » GR < 50 ) >GI (50> » OMEGA (999) *Y(10) »X(10) » 

1 NUMC(20) rKONTUO) *G(20f99) »G0R(999) »G0I(999) #NRATOR(10»20) * 

2 NDENOM ( 10 » 20) t GCOMR (5r 999) t GCOMI (5» 999) #pFXl(5»50) » 

3 PFY 1 ( 5 r 50 ) » KPTS(l) r XCOMN ( 10 » 50 ) » YCOMN ( 10 » 50 ) »KNR ( 1 ) #KDR( 1) 
DOUBLE PRECISION G 

IMPLICIT REaL*8 ( A-F » P-W ) 

REAL*4 X»Y» XV »CNR*CNI*XCOMN» YCOMN rCDR»CDI 
DIMENSION PPT (4) » LSN(4) 

COMMON TYPE ( 50 ) 

INTEGER TYPE 

COMPLEX P»PPT 

DO 140 U = 1 * NFREQ 

IF (TYPE (U) .EQ. ’G' )P=-PPT(1) 

IF ( TYPE (U) .EQ.’P* )P=-PPT(2) 

IF ( TYPE (U) .EQ. *S* )P=-PPT(3) 

IF ( TYPE (U) .EQ. 'A* )P=-PPT(4) 

IFITYPE(U) .EQ. *G» )SGN= LSN ( 1 ) 

IFITYPE(U) .EQ. *P» )S6N= LSN(2) 

IFCTYPE(U) .EQ. *S* )SGN= LSN(3) 

IF(TYPE(U) .EQ. * A* )SGN= LSN(4) 

KWHICH= KPTS(U) 

XV= OMEGA(KWHICH) * 6.2831853 
DO 130 L=1»KCHNL 
NCOMD= NUMC(L) 
knot=o 

LNOT=0 
I0P= KONT(L) 



DO 130 N=1»NC0MD 
IF(N.GT.l)l0P=2 
Nls NRATOR(l_»N) + 1 
Mis NDENOM(L»N) + 1 
DO 5 LP=1#N1 

5 X(LP)= XCOMN ( L » LP+KNOT ) 

DO 6 LPslrMl 

6 Y(LP)= YCOMN ( L » LP+LNOT ) 

K2=N1-1 

K3=M1-1 

CALL P0LFV(X»K2»XV»CNR»CNI) 

CALL P0LFV(Y»K3»XV»CDR#CDI) 

RD= CNR**2 ♦ CNI**2 

RR= (CDR*CNR+CDI*CNI)/RD 

RI=(-CDR*CNI+CDI*CNR)/RD 

GR(J)s GCOMR(L»KWHICH)*RR - GCOMI (LrKWHlCH) *RI 
GI(J)= 6C0MR(L#KWHICH)*RI + GCOMI (LrKWHICH) *RR 
As REAL(P)+GOR(KWHICH)-GCOMR(L#KWHICH) 

Bs AIMAG(P)+GOI(KWHICH) -GCOMI (LrKWHICH) 

FREQsl.o 

KSKlPsl 

DO 40 IslrNl 

KUL1=<-1)**( (I+D/2) 

KUL2s(-D**( (I+2)/2) 

FULlsKULl 

FUL2SKUL2 

IF(KSKIP-1>20>20»30 
20 C(I)s-GR(J)*FREG*FUL2 
D ( I ) s-GI ( J ) *FREQ*FULl 
KSKIP=2 
GO TO 40 

30 C(I)s-GI(J)*FREG*FUL 2 
D ( I ) s-GR ( J ) *FREG*FULl 
KSKlPsl 

40 FREQs FREQ*OMEGA(KwHlCH) *6.2831853 
FREQs l.o 
DO 50 IslrMl 
KMUL=(-1>**( (I+D/2) 

EMUL=KMUL 

E ( I ) s -FREQ * EMUL 

50 FREQs FREQ * OMEGA (KWHICH) *6.2831853 
FNAlsO.O 
FNA2S0.0 
DO 60 I=l»Nl 
FNAlsFNAl+C(I)*XCI) 

60 FNA2=FNA2+D ( I ) *X ( I ) 

FD2=0.0 

KIs 2 * ( (K3+1J/2) 

DO 70 I=2»Kl tZ 
70 FD2=FD2+E(I)*Y(I) 

FDlso.O 

KEs 2 * UK3 +2)/2) - 1 



DO 80 I = DKE»2 
80 FDl=FDl+Ed)*Yd) 

FN1= FNAl + FD1 * A - FD2 * B 
FN2= FNA2 + FD2 * A + FOl * B 
FD=FD1**2+FD2**2 
FN=FN1**2 +Fn2**2 

FYE= (FD *(A * FN1 + B * FN2 ) - FN * FDD/ FD**2 

FYO= (Fd*(”B * FN1 + A * FN2) - FN * FD2>/ F0**2 

FX1=FN1/FD 

FX2=FN2/FD 

PFX1(L»KN0T+1)= 0.0 

DO 90 I=1»KE»2 

PFY1(L»I+LN0T)= FYE * E(I) * SGN 
90 CONTINUE 

DO 100 I=2*KI»2 

PFY1(L»I+LN0T)= FYO * Ed) * SGN 
100 CONTINUE 

IF(I0P.EQ.2)PFYKL»LN0T+1)= 0.0 
DO 110 I=2»N1 

PFX1 (L* I+KNOT ) = (FXl*Cd) + FX2*Dd>) * SGN 
110 CONTINUE 

knot= knot + ni 

LNOT= LNOT + Ml 
130 CONTINUE 
KLAD=0 

DO 135 IX=1»KCHNL 
KNoT= KNRdX) 

DO 135 LXslrKNOT 
KLAD=KLAD+1 

135 G(jrKLAD)= PFXldX,|_X) 

DO 139 IX=1»KCHNL 
LNOT= KDRdX) 

DO 139 l_X=l * LNOT 
KLAD=KLAD+1 

139 G(O.KLAd)=PFY1(IX.LX) 

140 CONTINUE 
NPaRC=KLAD 

DO 150 U=1»NFREQ 
SUM=0.0 

DO 145 1=1 t NPARC 
145 SUM=SUM+G(J»I)*G(J»I> 

SUM= DSQRT(SUM) 

DO 149 I=1»NPARC 

149 G(j.I)= G (U» I ) /SUM 

150 CONTINUE 
RETURN 
END 


SUBROUTINE PHASEM (GR » GI r KPOINT » NM , FQMIN * FQMaX t KPTS » STBM . OMEGA) 
DIMENSION GR(1) . GI ( 1 ) » KPTS ( 1 ) . STBM ( 1 )' OMEGA (1 ) 
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C SUBPROGRAM FOR CALCULATING PHASE MARGINS 
C 

C DEFINITIONS OF I/O VARIABLES 
C 

C GR -ARRAY OF OPEN LOOP TRANSFER FUNCTION REAL PARTS 

C GI -ARRAY OF OPEN LOOP TRANSFER FUNCTION IMAGINARY PARTS 

C KPGINT-NO. OF POINTS 
C OMEGA -ARRAY OF FREQS. 

C NM -COUNTER 

C KPTS -FREQUENCY NOS. WHERE MARGINS OCCUR 

C STBM -STABILITY MARGINS OF MARGINS 

C FQMlN -LOWER FREQ. FOR MARGIN DETECTION 
C FQMAX - UPPER FREQ. FOR MARGIN DETECTION 
P=i.o 

DO 50 I=1»KP0INT 

SO= GR 1 1 ) **2 + GI ( I ) **2 

S2=SO-1.0 

IFII.EG.1)S1=S2 

IF (OMEGA ( I ) .LT.FQMIN) GO TO 40 

IF (OMEGA ( I ) .GT. FQMAX) RETURN 

IF(ABS(S2) .LT.1.0E-20)GO TO 30 

SGN=S2/ABS(S2> 

IF(S1*SGN.GT.O.O)GO TO 40 
30 11=1-1 

IF(A8S(S2) .LT.ABS(Sl) ) I 1=1 

nm=nm+i 

KP1S(NM) = U 

S3= (P+GRdl) >**2+Gl(Il)**2 
STBM(NM)= SQRT ( S3 ) 

40 S1=S2 
50 CONTINUE 
RETURN 
END 


SUBROUTINE GAINMG ( GR » GI » KPOINT ► NM » FQMIN » FQMaX * KPTS » STBM t OMEGA) 
DIMENSION GR(1) »GI(1) pKpTS(I) »STBM( 1) fOMEGA(l) 

C 

C SUBPROGRAM FOR CALCULATING GAIN MARGINS 
C 

C DEFINITIONS OF I/O VARIABLES 

c 

C GR -ARRAY OF OPEN LOOP TRANSFER FUNCTION REAL PARTS 

C GI -ARRAY OF OPEN LOOP TRANSFER FUNCTION IMAGINARY PARTS 

C KPOINT-NO. OF POINTS 
C OMEGA -ARRAY OF FREQS. 

C NM -COUNTER 

C KPTS -FREQUENCY NOS. WHERE MARGINS OCCUR 

C STBM -STABILITY MARGINS OF MARGINS 

C FOMIN -LOWER FREQ. FOR MARGIN DETECTION 
C FQMAX - UPPER FREQ. FOR MARGIN DETECTION 
P=1.0 
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DO 50 1 = 1 » KROINT 

S2=ei(I) 

IF(I,EQ.l)Sl=S2 
IF IOMEGA (I ) .LT.FQMIN) GO TO 40 
IF (OMEGA ( I ) *6T .FQM AX) RETURN 
IF(ABS(S2) .LT.1.0E-20)GO TO 30 
SGn=S2/ABS(S2> 

IF IS1*SGN.GT .0.0) GO TO 40 
30 CR= GR ( I ) 

IFlCR.6E.0.0)GO TO 40 
11 = 1-1 

IF(ABS(S2) .LT.ABS(Sl) ) 11=1 

nm=nm+i 

KPTS(NM)=I1 

S3= IP+GR(U) >**2+Gl (Il)**2 
STbM(NM)= SQRT(S3) 

40 S1=S2 

so continue 

RETURN 

END 


SUBROUTINE SRMINS (GR » GI > KPOINT »NM»P#N»FQMIN,KPTS»STBM»OMEGA) 
DIMENSION GR ( 1) »GI (1 ) rKpTS(l) »STBM(1) »OMEGA(l) 

SUBPROGRAM FOR DETERMINING THE MINMUNS OF THE OPEN LOOP FREQUENCY 
RESPONSE WITH RESPECT To THE -1 POINT WHEN GIVEN POINTS ON THE 
OPEN LOOP REQUENCY RESPONSE 

DESCRIPTION OF I/O VARIABLES 

GR - VECTOR OF REAL PARTS OF OPEN LOOP FREQUENCY RESPONSE 

GI - VECTOR OF IMAGINARY PARTS OF OPEN LOOP FREQ. RESPONSE 

KPolNT - NUMBER POINTS OF THE OPEN LOOP FREQ. RESPONSE GIVEN 
OMEGA - CORRESPONDING FREQUENCIES OF CHOSEN POINTS 
KPTS -FREQUENCY NOS. WHERE MARGINS OCCUR 
FQMIN -MINIMUM FRQ. CONSIDERED 
-P -POINT W.R.T. A MAX. OR MIN. IS DESIRED 

N -DETERMINES WHETHER A MAX, OR MIN. IS DETERMINED 

ASNl=0 . 0 
S1=0.Q 

DO 50 1=1 t KpOINT 
IF (OMEGA ( I ) .LE.FQMIN) GO TO 50 
S2= (P + GR ( I ) ) **2 + GI ( I ) **2 
ASN2=S2-S1 
5 CONTINUE 

IF l ASN2*N) 10»50’10 
10 IF ( ASNl*ASN2 )20»40»40 
20 IF ( ASNl*N ) 30 » 40 » 40 
30 NM=NM+1 
11= I - 1 
KPTS (NM) =11 
STbM(NM)= SQRT(Sl) 
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40 S1=S2 

ASNl= ASN2 
50 CONTINUE 
RETURN 
END 


SUBROUTINE DIRVEC (6 »NM»KPARC»DV» WEIGHT) 

directional vector program 
definitions of I/O variables 

g -matrix whose rows contain the gradient vectors of those 
stability margins only considered pertinent 
nm -number of stability margins considered pertinent 
weight-weighting factor vector 

DIMENSION G (20 » 99) » A(30>30)» WElGHT(l)» 

/ AI (30#30 ) t X ( 30 ) t DV ( 30 ) r Y(30) 

IMPLICIT REAL*8 ( A-H t O-Z ) 

DO 200 K=1'NM 
Y(K)= WEIGHT(K) 

DO 200 J=K»NM 
SUM= 0.0 

DO 150 I=1>KPARC 
150 SUM= SUM + G(Jrl) * G(K.I) 

A(jfK)= SUM 
A(K»U)= SUM 
200 CONTINUE 

IF (NM.GT .1 )GO TO 230 
AI(1»1)= 1«0/A(1»1) 

X(l)= WEIGHT(l) * AKl'l) 

GO TO 310 
230 CONTINUE 

CALL MATINV ( Af NM» AI » IER) 

IF(lER.EQ.0)GO TO 300 
WRITE (6»250 ) 

250 FORMAT (*0'»15X.*THE PARTIALS ARE NOT LINEARLY INDEPENDENT. 
/HE PROGRAM IS TERMINATED.') 

STOP 

300 CALL MATMUL(NM.AI»NM»Y»1»X) 

310 CONTINUE 

DO 450 I=1»KPARC 
SUM= 0.0 
DO 400 J=lrNM 

400 SUM= SUM + G( J» I ) * X(J) 

450 DV(I)= SUM 
690 RETURN 

end 


THUS T 


SUBROUTINE MATINV (2 . N t Y t IER ) 
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DO 50 I=1»KP0INT 
S2=GI(I> 

IF ( I ,EQ. 1 ) Sl=S2 
IF (OMEGA ( I ) .LT.FQMIN) GO TO 40 
IF (OMEGA ( I ) .GT .FQMAX) RETURN 
IF(ABS(S2) .LT.l.OE-20)GO TO 30 
SGn=S2/ABS(S2> 

IF (S1*SGN.GT *0.0) Go TO 40 
30 CR= GR ( I ) 

IF (CR.GE.O.O)GO TO 40 
11 = 1-1 

IFCABS(S2) .LT.ABS(Sl) ) I 1=1 

NM=NM+1 

KPTS ( NM) =11 

S3= (P+GR (II) ) **2+Gl ( II ) **2 
STtiM(NM) = SORT (S3) 

40 S1=S2 
50 CONTINUE 
RETURN 
END 


SUBROUTINE SRMINS(6R>GI»KP0INT#NM»P#N»FQMIN,KPTS*STBM»0MEGA) 
DIMENSION GR(1) rGId) »KPTS(1) fSTBM(l) »OMEGA(l) 

SUBPROGRAM FOR DETERMINING THE MINMUNS OF THE OPEN LOOP FREQUENCY 
RESPONSE WITH RESPECT To THE -1 POINT WHEN GIVEN POINTS ON THE 
OPEN LOOP REQUENCY RESPONSE 

DESCRIPTION OF I/O VARIABLES 

GR - VECTOR OF REAL PARTS OF OPEN LOOP FREQUENCY RESPONSE 

GI - VECTOR OF IMAGINARY PARTS OF OPEN LOOP FREQ. RESPONSE 

KPOINT - NUMBER POINTS OF THE OPEN LOOP FREQ. RESPONSE GIVEN 
OMEGA - CORRESPONDING FREQUENCIES OF CHOSEN POINTS 
KPTS -FREQUENCY NOS. WHERE MARGINS OCCUR 
FQMIN -MINIMUM FRQ. CONSIDERED 
-P -POINT W.R.T. A MAX. OR MIN. IS DESIRED 
N -DETERMINES WHETHER A MAX. OR MlN. IS DETERMINED 

ASNl=0.0 
S1=0.Q 

DO 50 1=1 .KPOINT 
IF ( OMEGA ( I ) .LE. FQMIN) GO TO 50 
S2= (P + GR(I))**2 + GI ( I ) **2 
ASN2=S2-S1 
5 CONTINUE 

IF ( ASN2*N ) 10 * 50 r 10 
10 IF(ASNl*ASN2)20»40f40 
20 IF ( ASNl*N ) 30 r 40 » 40 
30 NM=NM+1 
11= I - 1 
KP1S(NM)=I1 
STbM(NM)= SQRT(Sl) 
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MULTIPLIES (A) * (B) 

A IS AN NR X N 
B IS AN N X NC 
X IS AN NR X NC 

DO 4 1=1 rNR 

4 X ( I ) = 0.0 
DO 5 I=1»NR 
DO 5 K=1 *NC 
DO 5 J=lfN 

5 X(I)= X 1 1 ) + A ( I » J) * BCJ) 
RETURN 

END 


SUBROUTINE POLFV(Fw»K»X,FREAL*FIMAG) 

C PR06RAM FOR EVALUATING A POLYNOMIAL AT AN IMAGINARY FREQUENCY 
C 

c definitions OF I/O VARIABLES 
c 

C FW -VECTOR POLYNOMIAL COEFFICIENTS 
C K -ORDER OF POLYNOMIAL 

C X -FREQUENCY TO BE EVALUATED AT 
C FREAL -REAL PART OF Fw(JX) 

C FIMAG -IMAGINARY OF FW(UX) 

DIMENSION Fw(1) 

KEVEN=<K+2>/2 

K0DD=(K+l)/2 

Y=1 . 0 

freal=o.o 

DO 10 1=1 t KEVEN 

L=2*I-1 

FREAL= FREAL + FW(L)*Y 
10 Y=-Y*X*X 
FIMAG=0.0 

IF (KODD «EQ * 0 ) GO TO 30 
Y=X 

DO 20 1=1 »KODD 
L=2*I 

FIMAG= FIMAG + FW (L ) *Y 
20 Y=-Y*X*X 
30 RETURN 
END 


SUBROUTINE OTPT 1 ( STBM » OMEGA t KPT S > NM » XCOMN t YCOMN » KM IN » RQ » LOX » KCHNL » 
1 NUMC»NRAT0R»NDEN0M»PRX,PRY*I1>I2) 

DIMENSION STBMd) > OMEGA (1) rKPTS(l) > XCOMN ( 10 » 50 ) t YCOMN (10 #50) »RQ(1) 
1 » PRX ( 1 ) » PRY ( 1 ) # NUMC ( 1 ) t NRATOR ( 10 #20 ) • NDENOM ( 10 » 20 ) 

dimension type (so) 

COMMON TYPE 
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INTEGER TYPE 
C 

C SUBPROGRAM FOR THE OUTPUT OF INFORMATION CALCULATED 
C 

I0P=1 

WR1TE(6»10) LOX 

10 FORMAT (’O’ »25X »’ ITERATION NO. *»I4) 

DO 110 I=1»NM 
KWH= KPTS(I) 

FRLQ= OMEGA (KWH) 

IF ( I .EQ.KMlN+1 ) 60 TO 50 
IF ( I .EQ.l ) GO TO 70 
eo TO 90 
50 WRITE(6.60) 

60 FORMAT ( *0' »25X. ’ATTENUATED FREQUENCY INFORMATION’//) 

GO TO 90 
70 WRITE ( 6 r 60 ) 

80 FORMAT(’0» *25X» ’RELATIVE STABILITY INFORMATION’//) 

90 CONTINUE 

WRITE (6. 100 ) I. STBM(I)» FREQ. RQ(I)» TYPEd) 

ioo format ( * *»2X» ’margin radius no, * » i2 » *=• .fio.s.sx. *frequency=* . 

1 Flo. 5» IX, ’HZ’ ,5X. ’DESIRED MARGINz’ .FlO. 5, 5X» ’MARGIN TYPE=’»1X» 

2 Al) 

110 CONTINUE 

DO 104 I=1»KCHNL 
WRlTE(fa,102) I 

102 FORMAT(’0’»25X,* CHANNEL N0.’»I2.' COMPENSATORS’) 

L=NUMC ( I ) 

LX=1 

LY=1 

KX=0 

KY=0 

DO 104 IX=1.L 

KX=KX + NRATOR ( I . IX) + 1 

KY=KY + NUENOM(I.IX) + 1 

WRITE (6, 106) ( XCOMN ( I »N) »N=LX,KX) 

WRITE (6, 107) (YCOMN(I.N) »N=LY.KY) 

106 FORMAT ( ’ 0 ’ . 10X . ’NUMERATOR * »8F10 .5) 

107 FORMAT ( ’O’ »8X. ’DENOMINATOR* .8F10.5) 

LX=KX+1 

LY=KY+1 
104 CONTINUE 

WRITE (6. 130 ) (PRX(I) .I=I0P»I1) 

WRITE(6»120) (PRY(I) »I=I0PfI2) 

120 FORMAT ( ’0* . 'PARTIAlS WITH RESPECT TO Y’»8El0.3) 

130 FORMAT ( *0’ » ’PARTIALS WITH RESPECT TO X’.8El0.3) 

WRITE (6 » 160 ) 

160 FORMAT (’O’ ) 

RETURN 

END 



SUbROUTINE MULOUT (KCHNL. NUMC ► NRATOR »NDEN0 m> KNR .KDR» XCOMN» YCOMN) 
DIMENSION CON(30) . COM ( 30 ) . XCOF(30). XCOMN ( 10 . 50 ) . YCOMNtlO.50) 
1 NUMC (30) . NRATOR(10»20) . NDENOM<10.20) » KNR(20>. KDR(20) 

DO 80 I=1*KCHNL 
C0M1) = 1.0 

N=Q 

lax=i 

NAT= NUMC ( I ) 

WRITE(6.40) I 

40 format (* o*>25x .' channel no. '. i2»2x» ' compensator' ) 

DO 65 j=i»nat 

M=NRATOR(I.U) 

Ml— M + 1 
LAY= lax + M 
KL=0 

DO 62 K=LAX,LAY 
KL=KL+1 

62 COM(KL)= XCOMN(I.K) 

LAX= LAY + 1 

CALL POLMUL ( CON . COM . N . M » XCOF ) 

n=n+m 

N1=N+1 

DO 64 K=1.N1 

64 CON (K ) - XCOF(K) 

65 CONTINUE 
WRITE (6.67) 

67 format( 'O' .25X. 'Numerator coefficients*) 

WRITE ( 6 . 69) (CON( J) ,J=1.N1) 

69 FORMAT ( ' 0 ' » 2X » 7E15 . 5 ) 

CON(1)=1.0 

N=U 

LAX=1 

DO 75 U=1.NAT 
M= NDENOM(I.J) 

Ml= M+l 
LAY= LAX+M 
KL=0 

DO 72 K=LAX,LAY 
KL=KL+1 

72 C.OM(KL) = YCOMN(I.K) 

LAX= LAY + 1 

CALL POLMUL (CON. COM. N.M. XCOF) 

N=|\+m 

ni=n+i 

DO 74 K=l.Nl 

74 CON(K)= XCOF (K ) 

75 CONTINUE 
WRITE(6.77) 

77 FORMAT( *0* »25X. 'DENOMINATOR COEFICIENTS* ) 

WR 1TE (6 . 69) (CON(J) . J=1.N1) 

80 CONTINUE 
RETURN 
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END 


SUBROUTINE POLMUL (CON. COM.N.M. XCOF) 

DIMENSION CON(l). C0M(1). XCOF(l) » CONA(50) » COMRA<50) 

FOR DOUBLE PRECISION REMOVE C FROM FIRST COLUMN OF NEXT CARD. 
DOUBLE PRECISION CON. COM. XCOF. CONA. COMRA 

THE VECTOR CON IS A VECTOR OF THE COEFFICIENT OF A POLYNOMIAL 
OF ORDER N. 

THE VECTOR COM IS A VECTOR OF THE COEFFICIENTS OF A POLYNOMIAL OF 
ORDER M. 

THE VECTOR XCOF IS A VECTOR OF THE COEFFICIENTS OF THE PRODUCT OF 
A POLYNOMIAL OF ORDER N AND A POLYNOMIAL OF ORDER M. THE 
POLYNOMIAL OF WHICH THE COEFFICIENTS ARE THE VECTOR XCOF HAS AN 
ORDER OF M + N. 

DO 1 I=1»M 

1 CONAlI)=0.0 
NX=N+1 

DO 2 1=1. NX 

LX=M+I 

2 CONA(LX)=CON(I> 

MX=M+1 

DO 3 I=1»MX 
MY=M+2-I 

3 COMRA(I)=COM(MY) 

DO 4 1=1 » N 

NX=M+1+I 

4 COMRA(NX)=0.0 
KY=M+N+1 
KX=KY 

DO 7 K=1.KY 
XCOF(K)=0.0 
DO 5 L=1.KX 

5 XCOF (K)= CONA(L) * COMRA (L) +XCOF <K ) 

KX=KX-1 

DO 6 U=1 »KX 

6 CONA(U)=CONA(J+1) 

7 CONTINUE 
RETURN 
END 
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